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Preface 


The  results  of  comparing  real  and  synthetic  texture  during  the  recent 
Smart  Weapons  Operability  Enhancement  Joint  Test  and  Evaluation 
suggested  that  a  better  way  of  generating  texture  for  placement  in 
synthetic  scenes  be  developed.  In  this  study,  we  investigate  three 
accepted  measures  of  texture;  correlation  length,  fractal  dimensions, 
and  parameters  derived  from  the  grey-level  co-occurrence  matrix  as 
candidates  for  a  different  texture  generation  algorithm.  The  output  of 
these  measures  (or  metrics)  are  used  as  input  to  three  texture 
generating  algorithms.  A  comparison  is  then  made  between  the  input 
and  output  textures.  The  results  of  this  study  were  not  definitive  and 
future  work  is  indicated.  This  study  was  supported  by  the  Smart 
Weapons  Operability  Enhancement  Joint  Test  and  Evaluation  Program 
Office,  Dr.  J.P.  Welsh,  director. 
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Executive  Summary 


In  the  process  of  creating  synthetic  scenes  for  use  in 
simulations/visualizations,  texture  is  created  as  a  surrogate  for  high  spatial 
definition.  For  example,  measuring  the  location  and  characteristic  of  every 
blade  of  grass  in  a  lawn,  then  simulating  a  scene  of  it  would  be  excruciatingly 
laborious.  Various  techniques  have  been  devised  to  place  the  required  details 
in  the  scene  through  the  use  of  texturing.  Experience  gained  during  the  recent 
Smart  Weapons  Operability  Enhancement  Joint  Test  and  Evaluation  (SWOE 
JT«&E)  has  shown  the  need  for  higher  fidelity  texturing  algorithms,  and  a 
better  parameterization  of  those  that  are  in  use.  This  study,  analyzes  four 
aspects  of  the  problem:  texture  metrics,  texture  creation  algorithms,  texture 
extraction,  and  texture  insertion. 

The  overall  idea  is  to  see  if  a  textural  property  can  be  measured  with  a  metric, 
and  use  it  as  a  seed  for  the  creation  of  texture  with  that  same  textural 
property.  Textural  metrics  can  be  parameterized,  resulting  in  texture  created 
for  insertion  into  a  scene.  Texture  extraction  is  the  problem  of  determining 
metrics  capable  of  capturing  the  textural  aspects  of  a  scene  element.  Texture 
insertion  is  just  the  reverse  problem;  mapping  a  texture  into  a  scene  element. 
Furthermore,  if  the  texture  metrics  are  parameterized,  the  task  of  synthetic 
scene  generation  is  simplified.  For  example,  if  a  particular  metric  describes 
the  texture  of  a  certain  scene  element  and  that  metric  varies  only  with  time  of 
day,  then  that  parameter  (time),  can  be  used  to  determine  what  the  texture 
should  measure  at  that  time.  The  question  becomes,  can  we  generate  a  texture 
that  measures  properly  according  to  some  metric?  The  major  purpose  of  this 
study  is  to  determine  the  algorithms  that  use  certain  texture  metrics  as  seeds, 
and  to  generate  texture  maps.  These  maps  are  studied  to  determine  if  they 
yielded  correct  measures  of  their  texture.  The  three  texture  metrics  used  in 
this  study  are  those  used  in  the  SWOE  JT&E  analyses:  correlation  length, 
metrics  derived  fi-om  the  grey  level  co-occurrence  matrix,  and  fractal 
dimension.  Reviewers  of  this  report  noted  that  some  of  our  fractal 
dimensions  are  not  of  the  expected  value.  The  textured  metrics  used  are 
limited  in  their  use  and  application.  We  did  not  attempt  to  study  this  situation 
in  detail,  although  it  is  addressed  in  the  conclusions  of  the  report.  The  texture 
generation  algorithms  used  in  this  study  are  based  on  the  same  three  metrics. 
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1.  Introduction 


Haralick  and  Shapiro  define  texture  as  the  property  of  an  image  where  small 
regions  of  the  image  experience  “wide  variation  of  tonal  features.”  [1]  This 
is  opposed  to  the  grey-tone  property  where  small  regions  experience  “little 
variation  of  discrete  tonal  features.”  The  small  regions  make  up  the  set  of 
feature  and  background  fragments  of  an  image,  which  can  be  resolved  as 
unique  portions  of  the  image.  The  textural  properties  are  below  the  resolution 
threshold  for  definition  of  features.  A  grey-tone  feature  is  a  region  of 
contiguous  pixels  that  are  all  at  or  near  the  same  grey  level:  all  white,  all 
black,  or  somewhere  in  between.  A  textural  feature  is  a  region  of  contiguous 
pixels  whose  grey  level  varies  from  pixel  to  pixel,  or  from  a  couple  of  nearby 
pixels  to  another  couple  of  nearby  pixels.  The  variation  takes  on  a  continuum 
of  characteristics  from  a  variety  of  highly-ordered,  or  structured,  patterns  to  a 
variety  of  random  patterns.  The  textural  features  can  be  composed  of  large, 
grey-tone  features;  small,  grey-tone  features  or  a  combination.  The  measures, 
or  metrics,  are  the  tools  used  to  describe  texture  and  are  based  on  a  variety  of 
mathematical  operators  that  are  applied  to  the  2-D  array  of  pixels.  The  real- 
world  properties  that  cause  texture  to  exist  are  relief,  and  optical  and  thermal 
characteristics  (such  as  emissivity,  thermal  conductivity  of  the  various  scene 
elements),  which  exist  on  a  small  scale  relative  to  the  rest  of  the  elements 
(features)  in  the  scene. 

This  paper  addresses  how  the  fine  detail,  which  is  present  in  the  real  world, 
be  represented  in  synthetic  imagery  without  resorting  to  the  laborious  and 
time-consuming  efforts  of  including  the  details  in  a  modeling  and  simulation 
effort. 

Figure  1  provides  an  example  of  high-resolution  texture  map.  This  image 
was  obtained  during  the  SWOE  JT&E  [2]  by  the  Waterways  Experiment 
Station  (WES)  of  the  U.S.  Army  Corps  of  Engineers.  It  is  a  high  spatial 
resolution  image  of  a  small  region  of  natural  terrain  located  at  the  Yuma  I 
field  test  site  (one  of  the  locations  used  during  the  SWOE  JT&E  series  of 
field  tests).  The  data  was  obtained  with  a  far  infrared  (IR)  (8  to  12  |im), 
thermal  imager  configured  with  a  narrow  field-of-view  lens.  The  texture 
images  provide  the  raw  material  for  the  SWOE  JT&E  synthetic  texture 
generation  approach  to  the  problem  of  creating  fine-detail  synthetic  texture 
without  time-consuming,  highly-detailed  modeling. 
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Figure  1.  An  example  of  a  texture  image,  obtained  during  the 
SWOE  JT&E  to  define  a  texture  field  for  input  to  the  SWOE  image 
generating  process. 


Textural  measures  used  to  compare  the  real  and  synthetic  imagery  indicate 
that  the  procedures  for  introducing  texture  into  a  synthetic  image  are  flawed, 
in  some  cases  as  indicated  by  the  results  of  the  analyses  of  the  SWOE  JT&E 
synthetic  imagery.  [3]  Figure  2  illustrates  this  problem.  The  figure  consists 
of  three  panels:  a  real  image,  a  synthetic  image,  and  a  homogeneous  region 
map.  The  variation  in  grey  level  within  homogeneous  regions  (defined 
approximately  by  the  panel  on  the  right)  in  the  synthetic  image  is  generally 
lacking.  There  are  artifacts  in  the  various  panels  that  should  be  ignored  for 
the  purposes  of  the  discussion  here,  such  as  the  relative  contrast  between  the 
trees  and  terrain  for  real  and  synthetic  images,  and  the  border  around  the  real 
image. 
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Figure  2.  Three  panels  demonstrate  the  presence  of  texture  in  the  real  scene. 


Developing  texture-generating  algorithms  and  comparing  the  resulting 
synthetic  textures  among  themselves,  and  with  real  texture  is  the  focus  of  this 
paper.  The  textural  metrics  and  the  texture  generating  algorithms  used  in  this 
investigation  are  briefly  discussed.  A  comparison  is  made  between  the  output 
of  the  algorithms  and  the  real  textures. 


2.  Texture  Metrics 


Texture  metrics,  or  measures,  are  used  to  quantify  a  textural  property  of  the 
region  of  interest  so  that  various  regions  may  be  compared  to  see  if  they 
contain  similar  or  different  textures.  For  example,  in  the  SWOE  JT&E  Final 
Report  [2]  textural  metrics  were  used  to  determine  if  the  same  portions  of  the 
real  and  the  synthetic  image  were  actually  the  same.  The  standard,  first-order 
metrics,  such  as  average  grey  level,  do  not  give  any  information  as  to  the 
textural  structure  of  the  image.  The  second  order  metrics,  such  as  correlation 
length  do  allow  comparison  of  textures  between  images.  These  metrics 
become  the  eyes  that  allow  descriptive  comments  such  as  herring-boned 
structure,  fine-grained,  coarse-grained,  random,  and  isotropic  to  be  made 
about  a  textural  feature.  These  metrics  can  also  be  used  to  segment  an  image 
into  regions  of  similar  texture  through  application  of  a  particular  metric  to 
small  regions  of  the  image.  These  metrics  are  rarely  used  to  actually  compare 
regions  for  the  determination  of  a  degree  of  similarity  because  there  are  no 
standards  for  these  quantities  to  rank  the  texture  feature  under  study.  For 
example,  a  particular  feature  is  0.9  out  of  1 .0  of  having  the  particular  property 
being  measured. 

The  three  measures  discussed  in  this  paper  are: 

Correlation  length:  The  lag  where  the  spatial  autocorrelation  declines 
to  1/e  of  its  maximum  value.  Lag,  in  this  instance,  is  the  distance 
measure  defined  by  the  sampling  interval  in  space,  such  as  a  pixel. 

A  derivative  of  the  grey-level  co-occurrence  matrix  (GLCM):  The 
GLCM  measures  “the  dependence  between  pairs  of  grey  levels  arising 
from  pixels  in  a  specified  spatial  relation.  "  [1]  For  example,  the 
number  of  times  a  given  grey-level  pair  occur  side-by-side,  relative  to 
one  another,  at  separations  of  two  pixels,  three  pixels,  and  so  on  is 
examined.  A  variety  of  metrics  (entropy,  contrast,  etc.)  based  on 
various  moments  of  the  GLCM  have  been  developed  that  provide 
varying  degrees  of  parameterization  of  the  texture. 

Fractal  dimension:  For  the  purposes  of  this  paper,  the  fractal 
dimension  is  defined  as  the  parameter  D  in  the  following  equation: 
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D  =  n-\-\  —  H 


(1) 


with  «  =  2,  the  dimension  of  the  process  (a  line  has  a  value  of  1,  a  two- 
dimensional  process,  2,  and  a  three-dimensional  process,  3),  and 

2H  =  p-n.  (2) 

Under  the  assumption  that  the  image  grey-level  variation  can  be  modeled  as  a 
“fractional  Brownian  ‘noise’”  (0  <  /f  <  1),  the  value  of  p is  determined  from 
the  slope  of  the  power  spectral  density  function  of  the  image.  The  fractal 
dimension  then  measures  the  degree  of  roughness  of  the  image.  [4]  [5] 

These  metrics  are  used  to  compare  various  synthetic  textures.  These  metrics 
can  also  be  used  as  “seeds”  for  algorithms  to  generate  texture. 

There  is  some  question  as  to  some  of  the  results  of  applying  these  three 
metrics  to  a  particular  sample  of  data.  In  particular,  some  of  the  fractal 
dimensions  we  arrive  at  are  less  than  2.0.  According  to  the  literature,  this 
would  indicate  a  line  and  not  a  surface.  All  the  metrics  we  use,  and  many 
others  as  well,  do  not  behave  predictably  all  of  the  time  for  a  variety  of 
reasons.  In  the  case  of  the  GLCM,  the  way  an  image  is  quantized  (eight  bits 
versus  four  bits,  or  even  two  bits)  results  in  differing  GLCMs.  There  is  no 
guidance  for  choosing  one  over  the  other.  There  are  several  ways  to  define 
correlation  length.  Hilgers  et  al.  compare  different  measurements  of 
correlation  length  and  assess  the  stability  and  precision  of  their  results.  [6] 
They  conclude  that  their  results  depend  on  how  the  correlation  length  is 
computed,  which  depends  on  how  the  image  is  sampled,  and  the  finite  size  of 
the  sample. 

As  we  compute  it,  the  fractal  dimension  depends  on  the  measurement  of  the 
slope  of  the  power  spectral  density  curve,  which  depends  on  the  spectral 
leakage  caused  by  the  windowing  fimction,  the  presence  of  trends,  and 
aliasing.  This  says  nothing  about  the  assumption  that  power-law  scaling 
exists.  Some  of  these  issues  have  been  addressed  by  Austin  et  al.  and  Fox.  [7] 
[8]  Fox  states  that  the  Fourier  analyses  has  a  strong  statistical  foundation 
behind  its  application  that  is  often  overlooked.  On  the  other  hand,  the  fractal 
analysis  does  not,  and  it  appears  as  though  this  is  still  true  today.  Fox 
executes  an  empirical  study  between  fractal  dimension  and  power-law 
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frequency  spectra  to  conclude  that  the  relationship  is  non-linear,  except  in  the 
region  near  /3=2,  which  is  also  the  region  where  the  slopes  cluster  due  to 
spectral  leakage.  [7]  This  does  not  necessarily  solve  the  problem  of 
questionable  numbers,  but  serves  to  show  that  much  remains  to  be 
accomplished  in  this  field  and  that  a  more  in-depth  analysis  should  be 
performed  along  the  lines  of  that  Fox  demonstrates. 
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3.  Texture  Algorithms 


A  texture  map,  generated  from  a  fractal  dimension  algorithm  is  expected  to 
produce  a  measure  that  is  the  same  value  as  the  fractal  dimension.  If  that 
measure  were  capable  of  capturing  all  aspects  of  the  real  texture  and,  at  the 
same  time,  the  texture-generating  algorithm  based  on  that  measure  was  a 
capable  algorithm,  the  other  textural  metrics  (correlation  length,  etc.)  would 
also  be  expected  to  provide  the  same  measure  for  the  real  and  synthetic 
textures.  If  some  aspect  of  this  expectation  is  not  met,  then  further  study  will 
be  suggested.  It  is  not  expected  that  any  one  of  these  algorithms  will  be 
successful  100  percent  of  the  time.  Therefore,  we  chose  three  measures  to 
use  in  three  different  algorithms  to  generate  texture. 

The  mid-point  displacement  algorithm  is  a  method  of  generating  a  two- 
dimensional  fractional  Brownian  motion  or  noise.  This  approximate  method  is 
described  in  The  Science  of  Fractal  Images.  [9]  The  mid-point  displacement 
algorithm  produces  an  image  that  only  approximates  a  fractional  Brownian 
"motion”.  Another  algorithm,  based  on  this  same  concept,  does  a  better  job. 
[10]  See  appendix  A. 

The  GLCM  method  for  building  the  textures  presented  here,  is  based  on  the 
algorithm  described  in  the  paper  by  G.  Lohmann.  [11]  In  the  algorithm,  a  set 
of  four  grey-level  co-occurrence  matrices  provide  feature  vectors.  Lohmann 
reports  that  the  four  primary  co-occurrence  matrices  (horizontal,  vertical,  and 
left-  and  right-diagonal)  contain  sufficient  information  to  synthesize  texture 
that  closely  resemble  textures  from  remotely-sensed  images  of  the  Landsat/TM 
and  ERS- 1/AMI  satellite-borne  sensors.  Our  idea  is  to  determine  how  well  this 
algorithm  performs  on  our  type  of  imagery. 

For  a  description  of  the  correlation  length  algorithm  used  in  this  study,  see 
appendix  G  where  the  basic  methodology  is  laid  out:  [12] 

Numerous  approaches  to  texturing  were  reviewed  and  a 
simplified  two-dimensional,  autoregressive  (AR)  model  was 
selected.  It  uses  correlation  length  in  vertical  and  horizontal 
directions  and  brighmess  mean  and  standard  deviation  as  input 
parameters  for  a  kernel  (process  of  order  [1,1])-  The  model  was 
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developed  from  an  AR  model  used  in  the  U.S.  Air  Force  Infrared 
Modeling  and  Analysis  (IRMA)  image  modeling  system. 

An  empirical  approach  was  taken  because  of  the  lack  of  general 
theory  or  models  on  thermal  IR  texture.  Using  this  approach 
homogeneous  surfaces  (grass,  bare  soil,  trees,  tree  lines,  etc.)  of 
interest  were  imaged  for  the  times  (or  under  similar  conditions) 
the  synthetic  scenes  were  to  be  generated.  Textural  features  of 
these  surfaces  were  measured  and  input  to  an  AR  texture 
generator  program,  which  generated  isotropic  Gaussian  texture 
maps.  These  maps  were  used  by  the  rendering  software  to  texture 
the  polygons. 

The  raw  output  of  the  basic  program,  the  isotropic  Gaussian  texture  maps,  is 
smoothed  to  reduce  vertical  and  horizontal  striation  artifacts  for  long 
correlation  lengths.  Both  raw  and  smoothed  maps  are  used  in  the  study  reported 
here. 

Figure  3  shows  the  target  texture  attempted  to  be  reproduced  by  the  various 
algorithms  discussed  in  this  section.  The  map  is  a  32  by  32  image  (segmented 
from  the  image/texture  map  of  figure  1),  which  has  texture  properties  that  are 
given  in  table  1 . 


Figure  3a.  Texture  that  is  the  target  for  the  various 
generating  algorithms. 


Figure  3b.  The  histogram  of  grey  levels  for  the  target 
texture  image  of  figure  3a. 


The  results  shown  in  table  1  are  derived  from  the  application  of  the 
corresponding  metric  to  the  image.  The  table  headings  are  acronyms  for  the 
metrics  discussed  briefly  in  section  2.  With  the  exception  of  the  mean  and 
variance,  the  metrics  measure  the  second-order  statistics  of  the  images.  All 
have  been  used  in  studies  of  the  textural  properties  of  images.  More  details  on 
the  form  and  application  of  these  metrics  can  be  found  in 
Bleiweiss  et  al.  [3] 

In  tables  1,  2  and  3  MEAN  is  the  average  grey  level  of  the  texture  map  (ranges 
from  0  to  128),  VAR  is  the  variance  of  the  grey  levels  in  the  texture  map, 
FRAC  D  is  the  fractal  dimension,  ACL  is  the  autocorrelation  length,  CONTR, 
CORR,  ENTRPY,  and  HOMOG,  are  the  GLCM  metrics  called  contrast, 
correlation,  entropy,  and  homogeneity,  respectively.  The  GLCM  metrics  are 
somewhat  self-explanatory.  More  detailed  explanations  of  these 
metrics  are  in  appendix  C. 

Table  1.  Texture  measure  values  for  the  target  texture 

Mean  VAR  FRACD  ACL  CONTR  CORR  Entropy  HOMOG 

57.247  338.116  1.965  3.537  222.219  0.667  4.644  0.014 


Mid-Point  Displacement  Algorithm 

The  mid-point  displacement  algorithm  is  a  method  of  generating  a  two- 
dimensional  fractional  Brownian  motion  or  noise.  This  approximate  method  is 
described  in  The  Science  of  Fractal  Images.  [9]  Quoting  from  this  reference, 
but  modifying  the  discussion  to  fit  the  two-dimensional  aspects  of  our  problem, 
a  two-dimensional  fractional  Brownian  motion  is  defined  as  a  two-dimensional 
process  (a  random  field)  with  the  following  properties: 

The  increments  X{s^,sf)  are  Gaussian  with  zero 

mean  and  the  variables  s  and  t  are  positions  in  the  process; 
i.e.,  t  is  the  x-y  coordinate  of  the  point  under  discussion  and  s  is 
another  point  at  another  x-y  coordinate  at  some  distance 
removed  from  t. 


The  variance  of  the  increments  X(r,,t2)- depends 
only  on  the  distance 


and  is  proportional  to  the  -  th  power  of  the  distance,  where 
the  parameter  H  again  satisfies  0<  H  <\.  Thus, 

,r,)  -  4*,  )  “  [Zk  -  s,)") 

or 


Var(2r(r,,r2)-^(‘^1’'^2 


(3) 


Figure  4  demonstrates  the  values  of  the  random  field  at  two  points.  They  are 
separated  by  distance  d  and  denoted  by  A  and  B.  In  this  notation,  a  fi'actional 
Brownian  motion  has  the  property  that  the  difference  between  A  and  B,  A-B,  is 
Gaussian  with  mean  0  and  variance .  That  is, 

Var(.4  -B)  =  e{(^  -  5)' )  =  .  (4) 

The  mid-point  displacement  algorithm  begins  by  finding  the  comer  values  A, 
B,  G,  and  D  of  a  two  dimensional  array  (figure  5)  and  then,  at  each  stage,  finds 
values  at  other  positions  relative  to  these  in  the  following  steps: 

1)  interior 

2)  edge 

3)  interior 

4)  interior  (symmetric). 

The  points  in  step  four  above  are  symmetric  about  the  diagonal  to  those  in  step 
three. 


21 


Figure  4.  A  picture  of  the  relative  positions,  A  and  B,  within  the  random 
field,  separated  by  distance  d. 


B* 

• 

nf 

• 

G* 

• 

£>• 

Figure  5.  Graphical  aid  to  assist  in  defining  the  detailed  steps  of  the  mid¬ 
point  displacement  algorithm. 

The  initial  step  of  this  algorithm  obtains  values  for  the  four  comers  A,  B,  G,  and 
D  by  independently  sampling  a  Gaussian  distribution  with  zero  mean  and 
variance  figure  5).  It  is  easy  to  show  that  equation  (4)  is  satisfied  for 

any  two  adjacent  values  on  an  edge.  For  example,  A  and  B: 

Var(A  -B)=  e{(^  -  Byj=  E{A^}^2E{A£}  +  }  (5) 

Using  the  independence  of  A  and  B  we  have, 

Var(^  -B)  =  e[a^]  +  e[b'^  }  =  cr  .  (6) 

This  verifies  equation  (4 )  for  the  pair  of  points  A  and  B.  The  same  result 
obtains  for  pairs  {B,D),  (A,G),  and  (G,D). 

Equation  (4)  is  not  satisfied  for  the  diagonal  pairs  (A,D)  and  (B,G).  From  the 
independence  of^  and  D  and  the  definition  of^  and  D  we  have.  Just  as  above. 
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Var(^  -  D)  =  . 


(7) 


But  if  equation  (4)  were  to  be  satisfied  by  the  pair  A  and  D,  which  is  separated 
by  a  distance  of  dyll,  we  would  have  the  following  different  value  for 
Var(^-D)  : 


/  \1H 

Var(^  -  Z))  =  2cr^^f/-\/2j  .  (8) 

The  mid-point  displacement  algorithm  uses  two  recursive  steps.  One  is  used  to 
find  the  value  of  an  interior  center  point  m  (figure  5).  The  other  is  used  to  find 
the  value  of  an  edge  point,  such  as  that  labeled  e.  The  labels  found  in  figure  5 
are  used  in  the  following  discussion  to  explain  the  respective  recursive  steps. 

In  figure  5,  the  value  m  at  the  center  of  a  square  with  comer  values  A,  B,  G,  and 
D  is  determined  by  the  following  formula: 

(^A  +  B  +  G  +  D)  .qx 

m  =  - - -  +  s  (9) 


where  ^  =  an  independent  Gaussian  variate  with  zero  mean  and  variance 


VarCf)  =  ct’ 


(10) 


In  this  step,  the  variance  of  the  independent  variate  decreases  by  a  factor  of 
(1/2)"  .  With  this  assignment  for  the  variance,  equation  (4)  will  not  be 
satisfied  for  the  pair  {A,m)  even  at  stage  one  of  the  algorithm. 


In  figure  5,  the  value  e  at  an  edge  is  found  by: 

{A-^B  +  m) 
e  = - ^ 


(11) 


where  7  =  an  independent  Gaussian  variate  with  zero  mean  and  variance 
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(12) 


Var(7)  =  o-'(^|j  . 

Equation  (4),  therefore,  will  not  be  satisfied  for  the  pairs  (e,A)  and  (e,B). 

The  process  for  the  first  two  stages  is  illustrated  by  figures  6  and  7.  Figure  6 
illustrates  the  points  determined  by  stage  one.  The  center  point  is  labeled  one 
and  the  edge  points,  labeled  two  through  five,  are  determined  in  this  stage.  No 
interior  points  are  determined  in  parts  three  and  four  of  stage  one. 


Figure  6.  Depiction  of  stage  one  of  the  mid-point  displacement  algorithm. 
The  values  at  locations  one  through  five  are  chosen. 

Figure  7  shows  the  result  of  stage  two.  In  part  one  of  stage  two  the  interior 
points  six,  seven,  eight,  and  nine  are  found  in  that  order.  In  part  two,  the  edges 
10  through  17  are  found.  In  part  three,  interior  points  18  and  19  are  found. 
Finally,  points  20  and  21  are  determined. 


Figure  7.  Depiction  of  stage  two  of  the  mid-point  displacement  algorithm. 
The  values  at  locations  six  through  21  are  chosen. 


Figure  8  shows  an  example  of  a  texture  map  produced  by  this  algorithm.  The 
mid-point  displacement  algorithm  produces  an  image  that  only  approximates  a 
jfractional  Brownian  motion.  As  the  algorithm  produces  finer  grids,  the 
variance  of  the  difference  of  two  nearby  points  roughly  gets  smaller  as 
prescribed  by  equation  (4).  Appendix  A  describes  another  algorithm,  based  on 
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equation  (4).  Appendix  A  describes  another  algorithm,  based  on  this  concept 
[9]  which  does  a  better  job.  The  other  algorithm  was  not  discovered  until  this 
project  was  well  under  way;  otherwise,  it  would  have  been  afforded  a  more 
prominent  place  in  this  report. 


Figure  8a.  The  texture  image  output  by  the  mid-point  displacement 
algorithm. 


Figure  8b.  The  histogram  of  grey  levels  for  the  texture  image  of  figure 


8a. 
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3.2  Grey-Level  Co-Occurrence  Matrix  Algorithm 

The  Grey-Level  Co-Occurrence  Matrix  Algorithm  (GLCM)  method  for 
building  textures  presented  here  is  based  on  the  algorithm  described  in  the 
paper  by  G.  Lohmann.  [11]  The  GLCM  has  been  often  used  as  a  basis  for  the 
calculation  of  secondary  features  such  as  contrast,  correlation,  entropy,  and 
homogeneity,  which  are  used  to  categorize  and  segment  textural  features  in 
images.  In  this  algorithm,  a  set  of  four  GLCMs  is  used  to  provide  feature 
vectors.  Lohmann  reports  that  the  four  primary  GLCMs  (horizontal,  vertical, 
left-  and  right-diagonal)  contain  sufficient  information  to  synthesize  textures 
that  very  closely  resemble  textures  fi*om  remotely  sensed  images  of  the 
Landsat/TM  and  ERS- 1/AMI  sensors.  The  idea  is  to  determine  how  this 
algorithm  performs  on  our  type  of  imagery. 

Before  discussing  the  algorithm  we  present  some  basic  definitions.  Consider 
pairs  of  pixels  separated  by  distance  d  at  some  angle  (p.  Generally,  distances  of 
one  pixel  and  angles  of  0,  45,  90,  and  135°  are  used.  The 
{d=  1,  <!>=  0*’)-pairs  are  horizontally  adjacent,  the  (c?  =  1,  (j>=  OO^j-pairs  are 
vertically  adjacent,  the  (c/=  1,  ^  =  45°)-pairs  are  right-diagonal  neighbors,  and 
the  (d  =  1,  135°)-pairs  are  left-diagonal  neighbors.  If  n  denotes  the  number 

of  grey  levels  in  the  image,  then  the  (J,  ^-co-occurrence  matrix  Q  is  an  «  by  « 
matrix,  where  an  entry  (ij)  of  denotes  the  number  of  pairs  of  pixels 
separated  by  distance  d  at  angle  which  have  grey  values  i  andy. 

The  following  examples  of  0,  45,  90,  and  135°  co-occurrence  matrices 
(equations  (14)  through  (17))  for  the  image  called  “target”  below.  Equation 
(13)  illustrate  the  above  definitions: 

'o  1  o' 

target  image  =  2  12  (13) 

_0  2  1_ 

'O  1  1 

Co=  1  0  1 

0  2  0 


26 


(15) 


C45  = 


0  1  0 
1  0  0 
0  1  1 


^90  “ 


0  0  1 
0  1  1 
2  1  0 


,  and 


Q35  ~ 


0  0  0 
1  1  0 
0  1  1 


(16) 


(17) 


The  matrices  are  3  by  3  in  size  because  there  are  3  grey  levels  in  the  target 
image.  If  there  were  16  grey  levels  in  the  image,  then  the  GLCM  would  be  16 
by  16,  and  so  on. 

The  maimer  in  which  the  grey  level  co-occurrence  matrix  is  populated  follows. 
Referring  to  figure  9,  all  pairs  of  pixels  in  the  target  image  that  are  separated  by 
a  distance  of  1  pixel  at  an  angle  of  0°  are  determined  and  the  number  inserted 
into  the  0°  GLCM,  Co .  For  example,  the  0  1  situation  occurs  once; 

therefore,  a  1  is  placed  in  the  GLCM  at  the  position  (0,1)  as  shovra  in  panel 
(a)  of  figure  10,  while  the  2-^1  situation  occurs  twice  resulting  in  a  2  at 
position  (2,1)  as  shown  in  panel  (b)  of  figure  10.  Continuing  in  this  manner 
results  in  the  Co  given  by  equation  (14). 


0->1^0 

0  — ^  2  — ^  1 


Figure  9.  The  pairing  technique  used  to  populate  the  0°  GLCM  with  a 
distance  of  1  pixel  separating  pairs  in  the  horizontal  direction. 
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0  1  2 

0  1  2 

0 

1 

0 

1 

1 

2 

2 

2 

a 

b 

Figure  10.  Filling  the  0°  GLCM  of  the  target  image.  Panel  (a)  is  the  result  of 
counting  the  number  of  (0^1)  grey  levels  separated  by  a  distance  of  1  pixel  in 
the  horizontal  direction  (=1),  while  panel  (b)  is  populated  from  the  number  of 
(2-»l)  occurrences  (=2).  The  completed  matrix  is  given  in  equation  (14). 

1  0 

V'/ 

2  1 

Figure  11.  The  pixel  pairs  having  a  45°  relationship  and  whose  numbers 
are  used  to  populate  the  C45  matrix. 


Similarly,  for  the  45°  GLCM,  we  count  the  number  of  occurrences  of  pixel 
pairs  separated  by  one  space  in  the  45  direction  as  pictured  in  figure  1 1 . 
And,  as  with  the  0°  GLCM,  the  45°  and  remaining  GLCMs  are  given  in 
equations  (15)  through  (19). 

In  the  initial  step  of  the  algorithm,  a  random  image  is  generated  that  has  the 
same  grey  level  histogram  as  the  target  image  whose  texture  we  are  trying  to 
replicate.  The  histogram  is  calculated  from  the  row  sums  of  the  0, 45,  and  90 
degree  GLCMs: 

hist[z]  =  rs0[/]  -  rs45[/]  +  rs90[z]  +  tr[/] 


(18) 


where 


hist[/]  =  #  of  occurances  of  greylevel  i  in  the  "target"  image 
rsO[i]  =  sum  of  entries  of  row  /  of  the  0  degree  GLCM 
rs45[/]  =  sum  of  entries  of  row  i  of  the  45  degree  GLCM 
rs90[/]  =  sum  of  entries  of  row  i  of  the  90  degree  GLCM. 

The  vector  tr[/]is  the  top  right  comer  occupancy  vector.  That  is. 


1,  if  greylevel  i  occurs 
0,  otherwise 


at  the  bottom  right  comer  of "  target" 


The  values  of  these  vectors,  for  our  example,  are  given  by  the  following 
equations: 


rs0  = 


2 

2 

2 


9 


rs45  = 


1 

.1 

2 


rs90  = 


rsl35  = 


0 

2 

2 


,  and 


(19) 


(20) 


(21) 


(22) 


tr  = 


1 

0 

0 


(23) 
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Equation  (23)  can  be  verified  in  the  example  given  by  equation  (13).  To  find 
the  number  of  2's  in  the  target  image  we  calculate: 

hist[2]  =  rs0[2]  -  rs45[2]  +  rs90[2]  +  tr[2] 

=  2  -  2  +  3  +  0  (24) 

=  3. 

Note  that  the  indices  for  the  GLCM  matrices  begin  with  0.  The  formula  in  the 
computer  code  (appendix  C)  for  the  histogram  is  slightly  different  from 
equation  (23)  in  order  to  reflect  the  convention  of  placing  the  origin  of  an 
image  at  the  top  left  comer  rather  than  at  the  bottom  left  as  we  have  here. 

This  histogram  is  used  to  randomly  populate  a  matrix  to  produce  an  initial 
image.  Equation  (28)  is  an  initial  image  with  the  same  histogram  as  the 
example  given  by  equation  (13)  which  is  repeated  here  as  equation  (33). 
Equations  (29)  through  (32)  and  (34)  through  (37)  are  the  four  GLCMs 
associated  with  the  two  images.  The  sum  of  the  absolute  differences  of 
corresponding  positions  of  the  GLCMs  is  a  measure  of  the  closeness  of  the  sets 
of  the  GLCMs.  This  measure  of  differences,  known  as  the  Manhattan  metric, 
is  calculated  for  the  individual  GLCM  matrices  and  summed  to  19  in  equation 
(38).  It  should  be  noted  with  the  Manhattan  metric  there  are  a  variety  of 
distance  measures  used  to  describe  the  separation  of  two  vectors.  Three  of 
these  are  the  Euclidean  Distance,  the  Manhattan  (or  City  Block)  Distance, 
and  the  Square  Distance.  Mathematically,  these  are,  respectively: 


(25) 

(26) 

S-  MAX|A',-};|  . 

(27) 

In  the  case  of  the  Manhattan  Distance,  if  binary  vectors  replace  the  vectors 
and  i;,  then  the  distance  is  known  as  the  Hamming  Distance.  [13] 


initial  = 


(28) 


1  2  1 
1  0  0 
0  2  2 


Co  = 


1  0 
1  0 
0  1 


1 

1 

1 


Q5  = 


1 

0 

1 


1  0 
0  1 
0  0 


Qo  ~ 


0  2 
0  1 
2  0 


1 

0 

0 


0 

0 

1 


1  1 
0  0 
1  0 


target  = 


0  1 
2  1 
0  2 


0 

2 

1 


Co  = 


0  1 
1  0 
0  2 


1 

1 

0 


0  1 
1  0 
0  1 


0 

0 

1 


^90  ~ 


0  0 
0  1 
2  1 


1 

1 

0 


(29) 


(30) 


(31) 


(32) 


(33) 


(34) 


(35) 


(36) 


31 


(37) 


C  = 

'^135 


0  0  0 
1  1  0 
0  1  1 


=  4  +  6  +  3  +  6  (38) 

=  19. 

We  begin  a  sequence  of  pairwise  exchanges  in  the  initial  image  in  an  attempt 
to  arrive  at  an  image  that  is  closer  to  the  original  target  image  as  measured  by 
the  distance  metric.  If  the  distance  between  the  new  image  and  the  target 
image  is  less  than  the  distance  before  the  pairwise  exchange,  we  make  the 
exchange  that  results  in  a  new  image  in  the  iterative  process.  However,  if 
there  is  no  improvement,  we  may  make  the  swap  of  the  two  grey  level  values 
anyway.  We  can  call  this  a  "wildcard"  switch.  This  is  done  with  probability 
p  where: 


P  = 


1 

1  +  exp(^) 


(39) 


Here, 

A  =  the  current  error 
and 

T  =  a  temperature 

that  is  slowly  cooled  down.  This  method  is  a  way  of  simulating  an  annealing 
process  (see  appendix  B).  The  purpose  of  this  wildcard  switch  is  to  allow  exit 
from  a  “local  minimum.” 

The  algorithm  halts  if  the  number  of  successful  attempts  at  switching  falls 
below  some  predetermined  percentage  of  the  number  of  pixels  in  the  target 
image  lattice.  This  percentage  is  set  to  one  percent  or  less  of  M,  the  number 
of  pixels  in  the  target  image.  After  each  iteration  of  M  attempts  at  switching, 
the  procedure  either  halts  or  the  iteration  continues  with  a  new  cooler 
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annealing  temperature  T.  As  T  gets  smaller,  the  probability  p  also  gets 
smaller  and  eventually  there  are  fewer  and  fewer  "wildcard"  switches. 

When  two  grey  levels  are  swapped  in  an  image,  very  few  (eight,  at  most)  of 
the  entries  of  a  co-occurrence  matrix  are  affected.  The  entire  co-occurrence 
matrix  does  not  have  to  be  recomputed.  To  illustrate  the  effect  on  the  0 
GLCM  Co  by  way  of  example,  suppose  we  swap  the  two  grey  levels  that 
are  noted  by  an  asterisk  in  figure  12.  Then,  figure  13  shows  the  contribution 
of  these  grey  levels  to  the  GLCM,  and  figure  14  illustrates  the  changes  that 
occur  after  the  swap  has  occurred.  The  contributions  to  the  GLCM  by 
the  new  ordered  pairs  (1,0),  (0,1),  (1,2),  (2,0)  are  shown  in  the  right-hand 
side  of  figure  14.  The  new  Co  GLCM  is  obtained  by  subtracting  the  first 
change  matrix  fi-om  the  old  Co  GLCM  and  then  adding  the  second  change 
matrix,  as  illustrated  in  figure  15.  The  four  locations  in  the  new  Co  GLCM, 
which  have  changed,  are  flagged  by  an  asterisk.  The  two  locations  flagged 
by  a  double  asterisk  could  have  changed,  but  did  not  because  of 
cancellations.  In  the  most  general  case,  a  maximum  of  eight  locations  in  a 
GLCM  could  change  as  the  result  of  swapping  the  two  grey  levels. 

It  is  apparent  fi’om  this  example  that  it  is  importeint  to  keep  track  of  grey 
levels  that  are  at  the  tail  and  tip  of  the  arrows  in  figures  13  and  14.  The 
notation  that  is  used  in  the  computer  code  glcmtex.c  (appendix  G)  to  label 
these  points  is  shown  in  figure  16,  which  is  an  aimotated  version  of  figure  11. 
The  P-  denotes  the  “predecessor”  of  position  P,  the  grey  level  at  position  P-  is 
the  row  index  to  the  Cp  GLCM  entry  contributed  by  the  pair  I— >2.  Similarly, 
P+  indicates  the  successor  of  position  P.  The  grey  level  of  1  at  location  P+  is 
the  column  index  of  the  Q  GLCM  contribution  made  by  the  pair  2->l. 

In  terms  of  these  labels,  the  computer  code  glcmtex.c  considers  three  cases  in 
updating  the  GLCM  resulting  fi’om  a  swap: 

case  1  -  point  Q  coincides  with  point  P-; 

case  2  -  point  Q  coincides  with  point  P+; 

case  3  -  all  of  the  rest. 

These  labeling  conventions  are  used  in  the  fimctions  change  and  update_glcm 
in  glcmtex.c  and  are  explained  here  for  better  understanding  of  the  code.  The 
computer  code  also  checks  to  see  whether  points  P-,  P+,  Q-,  and  Q+  fall 
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outside  of  the  image.  For  example,  if  point  P  is  at  a  left  edge,  the  point  P-  does 
not  exist  when  considering  the  case  of  a  0°  GLCM. 

Figure  17  is  an  example  of  a  texture  map  produced  by  the  GLCM-based 
algorithm.  The  textural  properties  for  this  map  are  listed  in  table  2. 


initial  = 


12*1 
1  0*  0 
0  2  2 


Figure  12.  This  figure,  first  in  a  series  (12  through  15),  demonstrates  the 
affect  the  swapping  of  two  grey  levels  has  on  the  GLCM.  In  this  case,  the 
*-ed  values  are  the  two  positions  being  swapped. 


Figure  13.  The  second  figure,  of  the  series  (12  through  15),  indicates  the 
portion  of  the  0°  GLCM  that  is  affected  by  the  values  being  swapped,  and 
by  how  much  it  is  affected. 


34 


1  ^  0  ->  1 

1  2  0  yields: 

0  2  2 


+1 

+1  +1 
+1 


Figure  14.  The  third  figure,  of  the  series  (12  through  15 ),  indicates  the 
portion  of  the  0°  GLCM  affected  after  the  swap,  and  by  how  much. 
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Figure  15.  The  last  figure,  of  the  series  (12  through  15),  demonstrates  the 
technique  used  to  determine  the  new  GLCM  without  completely 
recalculating  the  whole  GLCM. 
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Figure  17a.  A  texture  map  created  through  the  use  of  the  GLCM-based 
technique. 


Figure  17b.  The  histogram  of  grey  levels  for  the  GLCM-based  texture 
map  of  figure  17a. 
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3.3 


Correlation  Length  Algorithm  (Current  Texture 
Generation  Process) 


The  two-dimensional  AR  model  used  in  this  study  is  found  in  the  report  by 
Sabol  and  Balick.  [14]  The  texture  generation  code  is  listed  there.  The 
PASCAL  code  was  translated  into  C  for  this  study,  and  is  included  in  appendix 
D. 

In  paragraph  6.4  of  appendix  G  of  the  SWOE  Final  Report,  [2]  the  basic 
methodology  is  laid  out: 

An  empirical  approach  was  taken  because  of  the  lack  of  general 
theory  or  models  on  thermal  IR  texture.  Using  this  approach 
"homogeneous"  surfaces  (grass,  bare  soil,  trees,  tree  lines,  etc.) 
of  interest  were  imaged  for  the  times  (or  under  similar 
conditions),  the  synthetic  scenes  were  generated.  Textural 
features  of  these  surfaces  were  measured  and  input  to  an  AR 
texture  generator  program,  which  generated  isotropic  Gaussian 
texture  maps.  These  maps  were  used  by  the  rendering  software 
to  texture  the  polygons. 

Numerous  approaches  to  texturing  were  reviewed  and  a 
simplified  two-dimensional  AR  model  was  selected.  It  uses 
correlation  length  in  vertical  and  horizontal  directions  and 
brightness  mean  and  standard  deviation  as  input  parameters  for 
a  kernel  (process  of  order  [1,1]).  The  model  was  developed 
from  an  AR  model  used  in  the  US  Air  Force  Infrared  Modeling 
and  Analysis  (IRMA)  image  modeling  system. 

Study  of  the  computer  code  reveals  that  texture  is  generated  one  pixel  at  a  time, 
using  the  following  equation: 


P[i,j]  =  al*  P[i,j  - 1]  +  a2*  P[/  - 1,;]  -  a\*  al*  P[i  - 1,7  - 1]  +  P[/,7]  (40) 
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where ; 


i  =  current  low  location 
j  =  current  column  location 
■P[^y]  =  pixel  at  location  [ij] 

R\i,j]  =  random  number  generated  at  location  [z,y] . 

The  random  number,  ,  is  constrained  such  that  following  conditions  are 
met: 


mean  =  m(\-al-a2  +  al*a2) 
variance  =  [l-a2*a2  — al*al(a2*al)^] 

where: 

m  =  mean  gray  level  of  the  image 
s  =  standard  deviation  of  gray  level  in  image 
a\  =exp(l//2) 
a2  =ext(l/v) 

h  =  horizontal  correlation  length  of  gray  level  variation  in  image 
V  =  vertical  correlation  length  of  gray  level  variation  in  image. 

Note  that  equation  (40)  is  not  the  same  as  the  corresponding  equation  on  page 
five  of  the  Sabol  and  Balick.  [14]  Equation  (40)  as  shown  here  reports  what 
is  in  the  PASCAL  code, 

Sabol  and  Balick  cite  that  several  modifications  were  made  to  the  original  AR 
model.  Modifications  are  listed,  with  our  comments: 

1)  Texture  used  in  this  application  is  isotropic  therefore  h  and  v 
are  set  to  the  same  value. 

2)  Mean  and  standard  deviation  are  fixed  at  128  and  32, 
respectively. 
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(In  our  study,  the  mean  and  standard  deviation  of  a  target- 
texture  image  are  used  as  input  parameters  to  the  computer 
program  to  fix  the  mean  and  standard  deviation  of  the 
generated  image.  The  model  automatically  produces  an  image 
with  the  specified  mean  and  standard  deviation.) 

3)  Because  the  AR  model  is  only  a  first-order  model  (based  on 
immediate  neighbors),  artifacts  (vertical  and  horizontal 
striations)  are  generated  for  long  correlation  lengths  (such  as, 
correlation  length  »  1).  A  moving-average  filter  (low-pass 
filter  )  is  applied  to  the  output  of  the  AR  module  to  reduce 
these  artifacts.  Only  the  diagonal  elements  are  included  in  the 
filters,  and  all  are  given  equal  weight.  The  size  of  the  filter  is 
set  to  the  correlation  length. 

(In  our  study  of  the  PASCAL  (procedure  moving-average), 
we  noted  that  the  low-pass  filtering  was  aecomplished  by  a 
convolution  filter  of  fixed  3  by  3  size.) 

4)  The  low-pass  filtering  described  above  compresses  the 
distribution  of  values  in  the  texture  map.  Larger  moving 
average  filters  result  in  greater  compression  of  the  distribution. 
The  desired  mean  of  128  and  standard  deviation  of  32  for  the 
distribution  is  restored  using  a  histogram  specification 
algorithm,  which  does  a  one-to-one  remapping  of  digital  values 
to  restore  the  desired  standard  deviation  to  32. 

(For  our  study,  the  remapping  for  the  filtered  image  was  not 
performed.) 

5)  Internally,  a  300  by  300  texture  array  is  generated,  but  only 
the  center  256  by  256  is  output.  This  eliminates  initialization 
effects  and  provides  a  large  enough  area  to  apply  the  moving 
average  filter. 

(For  our  study  a  32  by  32  pixel  image  was  generated.  No 
attempt  was  made  to  generate  a  larger  image  fi’om  which  to 
extract  a  32  by  32  subimage  was  made.) 


In  the  PASCAL  listing  for  the  function  RANI,  in  the  Sabol  and  Balick  report, 
the  constant  M2  is  defined  M2  =  124456.  [14]  This  varies  fi-om  Press  et  al. 


where  the  function  RANI  is  listed  with  a  PARAMETER  M2  =  134456.  [15] 
There  is  a  single  digit  discrepancy  between  these  two  numbers.  According  to 
Press,  the  choice  of  constants  used  to  implement  a  random  number  generator 
(of  which  this  parameter,  M2,  is  one)  is  critical  to  the  proper  operation  of  the 
generator.  However,  it  is  beyond  the  scope  of  this  report  to  delve  into  this 
aspect  further. 

Figure  18  shows  an  example  of  both  the  raw  and  smoothed  texture  map 
produced  by  the  algorithm.  The  textural  properties  of  these  maps  are  given  in 
table  2. 
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algorithm. 


4.  Comparisons  and  Discussion 

The  target  texture,  which  is  the  one  to  be  duplicated  by  the  various 
algorithms,  is  shown  in  figure  19  (a).  This  map  is  a  32  by  32  image 
(segmented  from  the  image/texture  map  of  figure  1)  that  has  texture 
properties  given  in  table  1 .  For  the  correlation  length  technique,  we  present  a 
typical  raw  image  (unsmoothed)  in  figure  19  (b).  In  figure  19  (c)  the  midpoint 
displacement  result  is  shown,  and  in  figure  19  (d),  the  grey  level  co-occurrence 
matrix  technique  output  is  demonstrated.  Table  2  contains  the  results  of 
applying  the  various  texture  metrics  to  these  texture  maps.  The  table  headings 
are  acronyms  for  the  metrics  discussed  in  section  2.  With  the  exception  of  the 
mean  and  variance,  these  metrics  measure  the  second  order  statistics  of  the 
images.  All  have  been  used  in  many  studies  of  the  textural  properties  of 
images.  The  target  texture  is  shown  in  figure  3,  the  texture  generated  by  the 
mid-point  replacement  algorithm  is  shown  in  figure  8,  figure  17  shows  the 
GLCM-generated  texture,  and  the  IRMA  (raw  and  smoothed)  is  shown  in 
figure  18.  Figure  19  repeats  many  of  these  maps  on  smaller  scale  for  side-by- 
side  comparison.  Reference  should  be  made  to  the  text  for  definition  of  the 
table  headings.  More  details  on  the  form  and  application  of  these  metrics  can 
be  found  in  Bleiweiss.  [3]  MEAN  is  the  average  grey  level  of  the  texture  map 
(ranges  from  0  to  128),  VAR  is  the  variance  of  the  grey  levels  in  the  texture 
map,  FRAC  D  is  the  fractal  dimension,  ACL  is  the  autocorrelation  length, 
CONTR,  CORR,  ENTRPY,  and  HOMOG,  are  GLCM  metrics  called  “contrast”, 
“correlation”,  “entropy”,  and  “homogeneity,”  respectively.  The  GLCM  metrics 
are  somewhat  self-explanatory  in  meaning. 


Table  2.  Texture  measure  values  for  generated  textures  shown  as 
examples  in  various  figures 


SOURCE 

MEAN 

VAR 

FRACD 

ACL 

CONTR 

CORR 

EIVTROPY 

HOMOG 

Target  Texture 

57.247 

338.116 

1.965 

3.537 

222.219 

0.667 

4.644 

0.014 

Mid-pt  replace 

56.740 

338.052 

1.774 

9.312 

22.554 

0.967 

6.641 

0.002 

GLCM 

57.247 

338.115 

1.804 

1.410 

376.695 

0.427 

4.776 

0.013 

IRMA 

58.790 

338.043 

2.117 

3.539 

164.708 

0.760 

7.180 

0.001 

IRMA(Smoothed) 

58.885 

211.740 

1.593 

4.314 

54.050 

0.873 

6.799 

0.001 
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Table  3.  Texture  measure  values  for  generated  textures 


Source 

MEAN 

VAR 

FRACD 

ACL 

CONTR 

CORK 

ENTRPY 

HOMOG 

Target  Texture 

57.247 

338.116 

1.965 

3.537 

222.219 

0.0667 

4.544 

0.014 

Mid-pt  replace 

average 

56.756 

336.775 

1.920 

11.116 

14.393 

0.978 

6.419 

0.002 

std.dev 

0.034 

4.268 

0.139 

8.274 

8.111 

0.012 

0.197 

0.00052 

GLCM 

average 

57.247 

338.116 

2.163 

1.452 

394.532 

0.399 

4.802 

0.013 

std.dev 

0.000 

0.000 

0.134 

0.402 

20.833 

0.032 

0.015 

0.00020 

IRMA 

average 

56.027 

327.618 

2.159 

3.055 

160.203 

0.748 

7.126 

0.001 

std.dev 

5.013 

67.931 

0.156 

0.811 

13.776 

0.045 

0.051 

0.00010 

IRMA 

average 

56.022 

213.424 

1.958 

3.813 

53.767 

0.864 

6.760 

0.001 

(smoothed) 

std.dev 

5.059 

63.775 

0.189 

1.211 

5.503 

0.031 

0.101 

0.00017 

Each  algorithm  was  run  101  times,  each  with  a  different  pick  from  the  random 
number  generator.  Each  algorithm  uses  a  random  number  generator;  however, 
the  specific  use  is  somewhat  different  from  algorithm  to  algorithm.  The  results 
of  these  calculations  are  shown  in  table  3,  which  are  the  same  as  table  2  except 
the  tabulated  values  applying  to  a  single  map  are  the  average  and  standard 
deviation  of  the  results  from  application  of  the  metrics  to  the  101  maps.  The 
sources  labeled  IRMA  and  IRMA  (SMOOTHED)  are  autocorrelation  length 
algorithm  oufr)ut.  The  variation  is  due  to  the  variation  introduced  by  the  use  of 
random  number  generators.  The  histograms  from  which  these  averages  and 
standard  deviations  were  determined  are  shown  in  appendix  D. 

Reference  to  table  3,  tells  the  story  about  the  relative  ability  of  the  various 
algorithms  to  replicate  the  target  texture  as  measured  by  the  various  metrics 
used  in  this  study: 

•  mid-point  replacement  algorithm  -  The  texture  map  created  by 
this  technique  is  similar  to  the  target  texture  as  measured  by  the 
fractal  dimension  (though  not  exact)  but  is  different,  as 
measured  by  toe  other  metrics. 

•  autocorrelation  length  (IRMA)  algorithm  -  This  texture  map 
(unsmootoed)  is  close  as  measured  by  toe  autocorrelation 
length,  but  is  quite  different  as  measured  by  toe  other  metrics 
(with  toe  possible  exception  of  toe  fractal  dimension  which  is 


similar  to  that  of  the  target).  The  smoothed  output  presents  an 
autocorrelation  length  much  closer  to  the  target.  The  fractal 
dimension  is  nearly  the  same  as  for  the  target,  in  fact  much 
closer  than  tiiat  of  the  mid-point  replacement  algorithm.  The 
GLCM  metrics  all  become  worse  under  the  action  of  the 
smoothing  function. 

•  GLCM  algorithm  -  The  texture  map  created  by  this  algorithm  is 
about  the  same  as  the  autocorrelation  length  algorithm  as 
measured  by  the  fractal  dimension  and,  except  for  the  entropy 
measure  (which  shows  good  agreement  with  the  target  entropy) 
all  of  the  other  GLCM  metrics  show  that  this  map  is  different 
from  the  target  map. 

It  would  appear,  then,  that  the  autocorrelation-based  technique  (at  least  the 
smoothed  version)  is  best  able  to  create  an  image  whose  autocorrelation  is 
nearly  the  same  for  target  and  synthetic,  while  the  other  techniques  seem  to 
not  be  able  to  replicate  themselves  very  well.  A  possible  exception  is  the 
IRMA  and  GLCM  ability  to  yield  a  fractal  dimension  close  to  that  of  the 
target.  It  remains  to  do  a  sensitivity  study  to  determine  the  discrimination 
ability  of  these  various  metrics  as  their  values  change  by  small  amounts.  For 
example,  is  there  much  difference  between  a  texture  whose  fractal  dimension 
is  2.75  and  one  whose  value  is  2.85,  and  so  on.  Visual  inspection  may  favor 
the  GLCM  algorithm,  though  this  may  be  due  to  the  level  of  pixellation  in  the 
two  maps  and  not  truly  comparable  texture. 

Why  then,  do  the  synthetic  SWOE  images  have  textural  metrics,  especially 
the  autocorrelation  length,  so  different  from  the  “real”  images?  [12]  More 
than  likely,  the  answer  lies  in  the  way  the  texture  was  measured  for  input  to 
the  SWOE  model,  the  way  in  which  the  texture  was  applied  to  the  synthetic 
images,  and  the  way  in  which  the  texture  was  measured  in  the  real  and 
synthetic  images.  For  example,  the  scene  that  was  imaged  for  textural 
characteristics  was  seen  at  a  different  perspective  view,  and  with  a  different 
spatial  resolution  than  that  of  the  real  images  and,  the  application  of  the 
synthetic  texture  is  made  to  the  scene  elements  as  seen  from  the  imager’s 
perspective  (but  with  a  texture  as  seen  from  a  different  perspective)  and  not  to 
the  scene  elements  before  they  are  rendered  into  an  image. 
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In  fact,  because  the  properties  that  create  texture  do  not  yield  imagery  that  is 
invariant  to  the  viewpoint,  it  is  not  at  all  clear  that  using  such  a  simple 
approach  would  work  at  all.  The  only  exception  is  a  simple  approach  where 
texture  is  extracted  from  the  real  images  and  then  pasted  into  the  synthetic 
images.  This  really  becomes  a  hybrid  technique,  not  a  true  synthetic 
technique.  For  an  invariant  approaeh,  the  texture  should  be  applied  as  a 
variation  of  small  scale  relief  and  optical  properties  (assuming  that  the 
thermal  properties  are  homogeneous  at  this  scale)  to  the  large  scale  relief 
Then,  when  the  seene  is  rendered  into  an  image,  the  texture  is  properly 
mapped. 

Future  work  may  also  look  toward  methods  to  parameterize  texture  based  on 
some  sort  of  environmental  parameter.  For  example,  one  approach  still  with 
many  limitations,  would  be  to  parameterize  the  measure  of  texture  for  a 
particular  imager  view  and  then  apply  the  correct  texture  based  on  the 
parameterization.  An  example  of  how  this  would  work  is  suggested  by  figure 
20,  a  plot  of  a  particular  texture  metric  with  time  of  day  for  the  Yuma  I  field 
exercise  of  the  SWOE  JT&E.  During  that  effort,  the  meteorological 
influences  were  quite  even  (day  to  day)  so  that  diurnal  trends  in  some  of  the 
textural  metrics  could  be  seen.  During  the  Grayling  field  exercises,  the 
meteorological  influences  were  so  variable  that  this  simple  correlation  was 
lost  to  ready  viewing;  though  it  may  possibly  be  extracted  through  some  sort 
of  multivariate  analysis. 


(c)  Mid-Point  Replacement  (d)  GLCM 


Figure  19a.  Four  textures  displayed  in  a  common  frame  for  comparison. 
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(e)  Target 


(f)  Correlation  Length 


(g)  Mid-Point  Replacement 


(h)  GLCM 


Figure  19  (b).  The  histograms  for  the  four  textures  of  panels  displayed  in  a  common 
frame,  for  a  better  visual  comparison  among  them. 


Figure  20.  Plot  of  the  entropy,  a  GLCM  metric,  versus  time  of  day  to  demonstrate 
the  diurnal  trend  experienced  by  a  metric  such  as  this,  which  could  be  used  in  a 
parameterization  of  a  seed  to  produce  synthetic  texture. 
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5.  Summary 


We  have  demonstrated  the  generation  of  texture  using  three  different 
algorithms  in  an  attempt  to  compare  their  various  abilities  to  replicate  a  given 
texture.  The  results  of  that  demonstration  are  mixed,  and  though  the 
autocorrelation  technique  was  certainly  no  worse  than  the  rest  of  the 
algorithms,  it  may  be  best  suited  for  the  task  at  hand.  We  feel  that  the 
synthetic  imagery  (as  judged  by  the  metrics)  generated  in  the  SWOE  process 
did  not  compare  more  favorably  with  real  image  as  to  the  textural  properties, 
because  of  the  way  the  texture  was  measured  and  then  applied  to  the  imagery. 
The  dependence  of  texture  on  view  angle  was  not  taken  into  account.  Also, 
the  texture,  as  seen  by  real  imagers  (for  the  case  studied  in  the  SWOE  effort) 
is  not  isotropic  nor  is  it  necessarily  homogeneous  over  those  scene  elements 
thought  to  be  homogeneous.  Definitive  answers  remain  until  a  more  in  depth 
study  can  be  performed.  The  first  step  has  been  taken;  the  generation  of  the 
software  and  methodology  to  perform  the  study. 
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6.  Future  Directions 


This  effort  has  focused  on  possible  quick  solutions  to  the  problem  of 
simulating  the  real  texture  of  a  scene,  in  lieu  of  modeling  the  scene  at  very 
high-spatial  resolution.  We  have  shown  that  there  does  not  appear  to  be  a 
quick  solution  to  the  problem  among  the  techniques  explored  here.  The 
reasons  for  this  include: 

•  texture  generating  algorithms  are  not  good  at  producing  a 
texture  with  the  same  measure  as  the  seed; 

•  texture  metrics  are  not  good  at  measuring  the  texture; 

•  some  of  the  texture  metrics  are  misused.  The  assumptions 
underlying  the  model  being  used  are  not  being  met. 

There  may  be  a  solution  of  the  type  we  are  looking  for  that  will  solve  the 
problem.  However,  our  research  provided  no  solutions.  Additional  work 
might  proceed  in  the  following  manner: 

•  further  inventory  of  available  texture  metrics  and  their 
appropriate  algorithms  for  use  in  generating  texture; 

•  grade  the  ability  of  these  metrics  to  discriminate  against  the 
type  of  textural  features  that  we  are  likely  to  encounter. 

Alternatives  range  from  the  very  simple  to  the  very  sophisticated.  For 
example,  choose  a  limited  variety  of  features  that  are  likely  to  generate 
different  textures,  model  them  in  detail,  and  then  view  them  under  a  limited 
number  of  conditions: 

•  sparse  grass  at  10  levels  of  sparseness; 

•  bare  dirt  at  different  levels  of  roughness; 

•  and  so  on. 

These  could  then  be  modeled  under  a  limited  number  of  seasons,  diurnal 
conditions,  and  look  angles  to  arrive  at  a  library  of  textural  features  that  could 
be  selected  to  fill  in  the  appropriate  polygons  of  the  scene. 


55 


To  help  in  this  modeling  effort,  we  need  to  acquire  a  library  of  features  that 
are  to  be  modeled.  This  could  include  synthetic  aperture  radar  (SAR)  data  to 
provide  the  relief  component  that  goes  into  the  model,  field  ground-truth 
measurements  to  ascertain  the  degree  of  homogeneity/inhomogeneity  in  the 
thermal  and  optical  properties  of  the  scene.  Hyperspectral  imaging  might  be 
useful  here. 

The  Rochester  Institute  of  Technology  [16]  as  well  as  many  other 
organizations,  such  as  Purdue  University  [17]  are  conducting  similar  efforts. 
The  Department  of  Defense  and  government  directors  of  research  such  as  the 
Army  Research  Office,  may  also  provide  pointers.  A  library  of  digital  texture 
maps  is  available  via  FTP,  over  the  internet  at 
Whitechapel. media. rait.edu/Vistex.  Other  studies  should  not  be 
overlooked,  as  they  may  be  of  some  help  in  ascertaining  the  utility  of  future 
investigations  by  providing  us  with  some  standard  textures. 
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Acronyms  and  Abbreviations 


AR 

autoregressive 

GLCM 

grey-level  co-occurrence  matrix 

IR 

infrared 

IRMA 

infrared  modeling  and  analysis 

SAR 

synthetic  aperture  radar 

SWOE  JT  &  E 

Smart  Weapons  Operability  Environment  Joint 

Test  and  Evaluation 

WES 

Waterways  Experiment  Station 
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Appendix  A 

Modified  Mid-Point  Displacement  Algorithm 


The  modified  mid-point  displacement  algorithm  is  a  method  of  generating  a 
two  dimensional  fi'actional  Brownian  "motion".  This  method  is  described  in 
Stoksik  et  al.  (1995).  Much  of  the  notation  and  wording  of  the  following 
description  comes  fi'om  this  paper. 


A  two  dimensional  fi-actional  Brownian  motion  (fBm)  is  a  random  process 
with  an  average  spectral  power  density  given  by 


s(Mi) 


(A-l) 


where  cr^  is  a  constant.  fBm  is  a  zero-mean  nonstationaiy  Gaussian  random 
process  B„{t)  indexed  by  the  parameter  The  symbol  t 

represents  the  vector  pair  (t,  ,^2 )  • 

The  e3q)onent  fi,  H,  and  the  firactal  dimension  5j  are  related  by  the 
equations: 

P  =  2H+2  (A-2) 

Sf=3-H  (A-3) 

The  statistical  behavior  of  fBm  can  also  be  described  by  its  covariance 
ejqjression, 

+lfr  -I' -  J|“)  (A-t) 

where,  for  instance, 

|f|  =  length  of  T 

and  s  ,  |?|,  and  |f-5|  are  defined  similarly. 

B„{i)  is  a  nonstationary  process  with  stationary  increments.  From  Equation 
(A-4)  we  can  derive  the  following  equation  asserting  the  finiteness  of  the 
structure  fimction  D{p) : 


D(p)  =  E^B„(t*p)-B„  (f )]’ )  =  2<7=  IpI’" 


(A-5) 


for  finite  p.  In  this  expression  it  is  evident  that  the  middle  expression 
involving  the  expectation  operator  E  is  independent  of  t.  The  purpose  of  the 
modified  random  displacement  algorithm  is  to  generate  a  sampled  grid  of 
points  where  the  statistical  interrelationship  of  the  points  is  given  by  Equation 
(A-5),  the  structure  function  condition.  Assuming  a  sampled  set  of  points  are 
available  which  satisfy  Equation  (A-5),  the  algorithm  generates  new  points 
between  the  existing  points.  These  new  points  are  generated  so  as  to  ensure 
that  they  have  the  correct  statistical  relationship  with  their  nearest  neighbors. 

The  originally  proposed  random  midpoint  displacement  algorithm  (RMDA) 
is  only  accurate  for  the  case  of  simple  Brownian  motion  {H  =  0.5).  For  other 
values  of  there  are  correlations  between  the  existing  random  variables  that 
must  be  considered  when  calculating  the  displacements  for  the  interpolated 
samples.  In  order  to  generate  accurate  fBm  with  the  correct  structure 
function,  the  algorithm  must  be  modified  to  account  for  the  correlations 
between  points  within  the  fflm. 

The  modified  mid-point  displacement  algorithm  begins  by  finding  the  comer 
values  of  a  two  dimensional  array  and  then  stage  by  stage  finds  the  values  at 
other  points  in  the  maimer  illustrated  by  Figures  A-1  and  A-2.  A  broad  view 
of  the  algorithm  is  gained  by  examining  the  organization  of  the  computer 
code  accfrac.c  (Appendix  I).  The  computer  code  shows  that  each  stage  of  the 
algorithm  determines  (1)  some  interior  points,  (2)  the  edge  points,  (3)  some 
more  interior  points,  and  then  (4)  interior  points  which  are  at  positions 
symmetric  to  the  diagonal  relative  to  those  determined  in  (3). 

The  modified  algorithm  starts  with  the  calculation  of  the  four  comer  samples 
denoted  by  A,  B,  G,  and  D  in  Figure  A-3.  If  these  four  comer  variates  satisfy 
Equation  (A-5)  we  have  the  following  six  equations  corresponding  to 
Equation  (A-5)  applied  to  the  six  possible  combinations: 
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E{(^-5)^}  =  2crV" 

(A -6) 

E{(^-G)^j  =  2o-V" 

(A -7) 

e{(5-I>)^}  =  2<tV" 

(A -8) 

E{(G-Z))^j  =  2cTV^" 

(A -9) 

E{(^-D)'}  =  2fT^(i/V2)'" 

(A -10) 

E^(B-Gf'^  =  2(T^[d42Y 

(A-11) 

The  last  two  equations  reflect  the  fact  that  the  distance  between  the  comer  sample 
locations  is  d^[2 . 

We  assume  that  the  values  A,  B,  G,  and  D  are  derived  from  the  following 
combination  of  the  auxiliary  variates,  S,  T,  U,  V,  W,  and  X: 


A  =  S  +  WI2 

(A -12) 

B  =  T+X/2 

(A -13) 

G  =  U-X/2 

(A -14) 

D  =  V-WI2 

(A -15) 

A  4 
2  1 
G  5 


B 

3 

D 


interior:  1 
edge:  2,3, 4,5 
interior:  none 
interior:  none 


Figure  A-1.  This  figure  illustrates  the  values  found  in  stage  1.  Only  the  center 
point  labeled  1  and  the  edge  points  2  through  5  are  determined  in  stage  1.  No 
interior  points  are  determined  in  parts  (3)  and  (4)  of  this  stage. 
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A 

12 

4 

16 

B 

10 

6 

20 

8 

11 

2 

18 

1 

19 

3 

14 

7 

21 

9 

15 

G 

13 

5 

17 

D 

interior:  6, 7,8, 9 
edge:  10-17 
interior:  18,19 
interior:  20,21 


Figure  A-2.  This  figure  shows  the  result  of  stage  2.  In  part  (1)  of  stage  2  the 
points  6,7,89  and  9  are  found  in  that  order.  In  (2),  the  edges  10  through  17.  In 
part  (3)  the  points  18  and  19  are  found.  Finally  points  20  and  21  are  determined. 


B* 

• 

• 

G* 

• 

D* 

Figure  A-3.  Graphical  aid  to  assist  in  defining  the  detailed  steps  of  the  Modified 
Mid-point  displacement  Algorithm. 

Here  S,  T,  U,  and  V  are  zero  mean  Gaussian  variates  with  the  same  variance  cr] 
yet  to  be  determined.  Both  W  and  X  are  also  zero  mean  Gaussian  variates  with 
the  different  variance  cr]  also  yet  to  be  determined.  The  variates  W  and  X 
force  correlations  between  the  comer  pairs  (A,D)  and  (B,G). 

Equations  (A-6)  through  (A-1 1)  force  conditions  on  the  two  variances  of  the 
six  auxiliary  variables  which  can  be  solved  to  yield  the  following  expressions 
for  the  values  of  a]  and  a]  : 

CTl=cr^(2-2'')rf“  (A-16) 

<T5=4<T'(2"-l)rf>"  (A -17) 
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The  modified  mid-point  displacement  algorithm  uses  two  recursive  steps. 
One  is  used  to  find  the  value  of  an  interior  center  point  m  as  in  Figure  A-3. 
The  other  is  used  to  find  the  value  of  an  edge  point  also  labeled  e  in  Figure  A- 
3.  The  labels  found  in  Figure  A-3  will  be  used  in  the  discussions  that  follow 
which  explain  the  respective  recursive  steps. 

In  Figure  A-3,  the  value  m  at  the  center  of  a  square  with  comer  values  A,  B, 
G,  and  D  is  determined  by  the  following  formula: 

(A  +  B  +  G  +  D^  -  a  1 0'! 

m  =  ^ - L  +  s  (A-18) 

4 


where  ^  is  an  independent  Gaussian  variate  with  zero  mean  and  variance 
Var(f).  The  value  of  Var(f)  is  obtained  as  follows.  From  Equations  (A- 12) 
through  (A-18),  we  can  determine 


(A-19) 


whereas  that  predicted  fi-om  Equation  (A-5),  the  stmcture  fimction  condition, 
is 


E{(A-r,y]=2.^[^y  . 

The  variance  of  the  displacement,  Var(f),  is,  therefore. 


(A-20) 


(A-21) 


This  value  for  Var(£)  will  also  force  Equation  (A-5)  to  hold  for  the  other  three 
pairs  (G,m),  and  (D,m).  By  comparison,  the  unmodified  RMDA  takes 
no  account  of  the  correlations  of  the  comer  samples  and  would  use  the 
variance  . 


Using  the  notation  of  Figure  A-3,  the  value  e  at  an  edge  is  found  by: 


(A-22) 


{A^B) 


where  77  is  an  independent  Gaussian  variate  with  zero  mean  and  variance 


(A-23) 


This  value  for  Var(;7)  is  found  through  evaluation  of  E{(A-e)^}  using  the 
definition  of  e  and  the  structure  function  condition.  By  way  of  contrast,  the 
RMDA  would  assign  the  following  value  for  the  variance  of  rj: 


Var(;;)  =  cT^[|j  =^^(2--)^-, 


(A.24) 


In  Stoksik,  et  al  (1995),  it  is  shown  that  this  method  yields  images  whose 
structure  fimctions  are  a  good  approximation  to  the  ideal  structure  functions. 
The  structure  functions  were  calculated  by  radial  averaging  of  the  two- 
dimensional  image  (surface)  and  then  plotting  the  resultant  one-dimensional 
curve  on  a  log-log  plot.  It  is  not  perfect  because  each  interpolated  point  is 
calculated  from  only  the  four  adjacent  points  surrounding  it  and  consequently 
the  structure  function  between  an  interpolated  point  and  a  more  distant  point 
is  not  necessarily  exact.  Figure  A-4  demonstrates  the  results  of  using  this 
algorithm  for  the  values  of  H  =  0.1  and  0.9,  respectively.  Figure  A-5  shows 
similar  results  for  the  original  mid-point  displacement  algorithm  as  discussed 
in  the  main  body  of  this  report. 
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Figure  A-6.  Histograms  of  grey  levels  for  the  images  of  Figures  A-4  and  A-5,  respectively.  The  upper  plots  are  for  the 
modified  mid-point  displacement  method  and  the  lower  panels  are  for  the  original  mid-point  displacement. 


Reference:  Appendix  A 


Stoksik,  M.A.,  R.G.  Lane,  and  D.T.  Nguyen,  “Practical  Synthesis  of  Accurate 
Fractal  Image,”  Graphical  Models  and  Image  Processing,  57  ,  3,  p  206-219, 
1995. 
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Appendix  B 


Simulated  Annealing:  A  Brief  Discussion 


The  primary  references  used  in  the  following  brief  discussion  are  Kirkpatrick, 
et  al.  (1983),  Kirkpatrick  (1984),  and  Metropolis  (1953).  The  McGraw-Hill 
Dictionary  of  Scientific  and  Technical  Terms  (Parker  1989)  defines 
annealing  in  the  following  way: 

“To  treat  a  metal,  alloy,  or  glass  with  heat  and  then  cool  to  remove 
internal  stresses  and  to  make  the  material  less  brittle.” 

Simulated  annealing  is  a  process  whereby  a  “system”  is  brought  to  a  “reduced 
energy  state”,  or  “temperature”  through  emulation  of  a  naturally  occurring 
process  -  this  is  accomplished  by  a  variety  of  algorithms.  The  goal,  in  our 
situation,  is  to  find  an  absolute  minimum  through  a  process  that  does  not  get 
led  astray  by  local  minima  nor  is  caught  for  good  in  a  local  minimum  --  we 
want  to  find  a  rearrangement  of  pixels  such  that  the  resulting  image  has  the 
same,  or  nearly  the  same,  GLCM  as  the  target  arrangement  of  pixels.  To 
accomplish  this  goal,  simulated  annealing  is  the  process  of  choice. 

In  metallurgy,  a  system  —  a  metal  —  is  heated  to  the  fluid  state,  and  the 
temperature  at  which  that  occurs  is  noted.  The  temperature  is  then  reduced, 
according  to  some  schedule  (the  “annealing  schedule”)  until  the  metal  has 
cooled  to  some  predetermined  state  --  usually  such  that  the  end  state  is  of 
uniform  condition;  i.e.,  no  cracks,  etc.  For  the  statistical  mechanics 
analogue,  a  system  of  particles  is  studied  at  “low  temperature”  which  state 
has  been  achieved  by  some  annealing  process.  What  happens  in  this  case,  is 
that  the  system  is  at  some  initial  configuration,  {x,},  with  some  potential 
energy  E  which  is  determined  by  the  Boltzmann  probability  factor: 

exp[-£({x,.})/A:r]  . 

For  our  purposes,  instead  of  the  energy  state,  we  use  a  “cost  factor’  defined  to 
be  A  and,  instead  of  just  decreasing  T  to  arrive  at  an  end  state,  we  use  what  is 
known  as  the  “Metropolis  Algorithm”  to  keep  the  iterative  process  moving. 
This  algorithm  simulates  the  behavior  of  a  many-body  system  at  some 
temperature.  Use  of  this  algorithm  then  “provides  a  natural  tool  for  bringing 
the  techniques  of  statistical  mechanics  to  bear  on  optimization.”  Still  quoting 
fi-om  Kirkpatrick,  et  al.  1984: 


Iterative  improvement,  commonly  applied  to  such  problems, 
is  much  like  the  microscopic  rearrangement  processes 
modeled  by  statistical  mechanics,  with  the  cost  function 
playing  the  role  of  energy.  However,  accepting  only 
rearrangements  that  lower  the  cost  function  of  the  system  is 
like  extremely  rapid  quenching  from  high  temperatures  to  T 
=  0,  so  it  should  not  be  surprising  that  resulting  solutions  are 
usually  metastable.  The  Metropolis  procedure  from 
statistical  mechanics  provides  a  generali2ation  of  iterative 
improvement  in  which  controlled  uphill  steps  can  also  be 
incorporated  in  the  search  for  a  better  solution. 

The  simulated  annealing  process  consists  of  first  ‘melting’ 
the  system  being  optimized  at  a  high  effective  temperature, 
then  lowering  the  temperature  by  slow  stages  xmtil  the 
system  ‘fi^ezes’  and  no  further  changes  occur.  At  each 
temperature  ,  the  simulation  must  proceed  long  enough  for 
the  system  to  reach  a  steady  state.  The  sequence  of 
temperatures  and  the  number  of  rearrangements  of  the  {x, } 
attempted  to  reach  equilibrium  at  each  temperature  can  be 
considered  an  annealing  schedule. 

Annealing,  as  implemented  by  the  Metropolis  procedure, 
differs  from  iterative  improvement  in  that  the  procedure  need 
not  get  stuck  since  transitions  out  of  a  local  optimum  are 
always  possible  at  nonzero  temperature. 

The  monograph  by  Beale  and  Jackson  (1991)  discusses  this  process  as 
applied  to  what  are  known  as  “Hopfield  Networks”  --  a  neural  network  which 
“...consists  of  a  number  of  nodes,  each  connected  to  every  other  node:  it  is  a 
fully-connected  network...”.  The  Hopfield  net  uses  simulated  annealing  to 
allow  exit  from  local  minima  —  large  changes  are  implemented  at  “high” 
temperatures  in  search  of  a  path  to  a  global  minimum.  The  temperature  is 
reduced  as  the  process  proceeds  which  allows  less  drastic  changes  to  be  used 
with  the  more  probable  result  that  the  process  ends  up  in  a  lower  energy  state. 
The  probability  distribution  which  controls  this  process  is  the  Boltzmann 
distribution. 
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Appendix  C 


Grey  Level  Co-Occurrence  Matrix:  A  Brief  Discussion 


Grey-level  co-occurrence  analysis  "...  characterizes  the  micro-texture  of  an 
image  region  by  measuring  the  dependence  between  pairs  of  grey  levels 
(our  grey  levels  are  in  units  of  apparent  temperature)  arising  from  pixels  in 
a  specified  spatial  relation"  [Haralick  and  Shapiro,  1991].  This  technique 
may  be  applied  to  either  the  whole  image  (global  feature  definition)  or  to  a 
ROI  (local  feature  definition).  Although  there  are  an  extremely  large 
number  of  parameters  to  be  derived  from  this  analysis,  it  can  be  limited  to  a 
few  that  have  been  shown  to  provide  good  discrimination  between  various 
textures  [Wahl,  1987;  Ballard  and  Brown,  1982;  Marceau  et  al.,1990]: 

1)  contrast  -  this  parameter  is  greatest  when  adjacent  pixels  are  very 
different  in  grey  level 

2)  homogeneity  -  this  is  greatest  when  most  of  the  co-occurrences 
are  for  the  same  two  grey  levels,  it  is  a  measure  of  how  consistently 
a  pattern  is  repeated  in  an  image 

3)  entropy  -  this  measures  the  "information  content  of  the  image;  it 
is  greatest  when  all  co-occurrence  possibilities  occur  in  equal 
proportion 

4)  correlation  ~  correlation  is  greatest  when  the  pixels  in  the  (i,j) 
pair  are  similar  in  grey  level  value  and  both  values  at  the  same  time 
are  either  above  or  below  the  average  values  for  their  respective 
positions  in  the  pair  (//,  and//^) 

Other  "texture"  analysis  techniques  are  also  available  [He  and  Wang,  1990; 
Therrien  et  al.,  1986]. 

The  grey-level  co-occurrence  matrix  (GLCM)  is  a  square  matrix  of  width  m 
based  on  the  number  of  grey  levels  in  the  image  under  analysis.  For 
example,  an  image  quantized  at  256  levels  causes  a  GLCM  of  256  x  256  to 
be  formed.  Because  this  creates  an  extremely  large  matrix,  the  number  of 
levels  is  generally  reduced  to,  generally,  8  or  16  [Wahl,  1987].  In  addition, 
as  will  be  seen  shortly  in  the  definition  of  how  the  GLCM  is  formed,  it  is 
possible  to  create  a  matrix  for  a  large  number  of  angles  or  orientations  as 
well;  this  too  is  usually  reduced  to  a  small  number;  0°,  45°,  90°,  and  135°. 


The  matrix  is  then  a  "distribution"  of  the  number  of  times  that  certain 
configurations  of  brightness  levels  occur.  This  can  best  be  explained  by 
example  [Haralick,  1974]: 

digital  image  (1  columns  and  k  rows): 

1=1234 

k=  1  0011 

2  0  0  1  1 

3  0  2  2  2 

4  2  2  3  3 

generalized  GLCM: 

grey  level 

0  1  2  3 

0  #(0,0)  #(0,1)  #(0,2)  #(0,3) 

1  #(1,0)  #(1,1)  #(1,2)  #(1,3) 

grey  level 

2  #(2,0)  #(2,1)  #(2,2)  #(2,3) 

3  #(3,0)  #(3,1)  #(3,2)  #(3,3) 


where  #(ij)  is  the  number  of  times  that  the  i  j  pair  occur  in  the  digital  image 
in  a  particular  orientation.  Specifically,  for  (t)=0°  (horizontal  orientation), 
there  are  four  times  that  a  zero  grey  level  value  occurs  left/right,  right/left  in 
the  digital  image  (in  some  instances,  only  the  left/right  or  only  the  right/left; 
i.e.,  uni-directional,  direction  is  chosen  which  would  yield  a  value  of  2  for 
this  case): 


79 


k,l=  1,1  is  next  to  1,2 
=  1,2  is  next  to  1,1 
=  2,1  is  next  to  2,2 
=  2,2  is  next  to  2,1. 

There  are  no  more  k,l  positions  where  a  zero  grey  level  occurs  and  right 
next  to  it,  another  zero  occurs.  This  process  continues  in  this  manner  to 
yield  the  0®  GLCM  : 

And,  we  can  continue  with  the  example  so  that  <|)  =  45°,  90°,  and  135°  are 
processed: 

Similar  matrices  can  be  obtained  for  distance  measures  greater  than  the 
value  of  1  which  we  have  used  —  we  could  have,  for  example,  determined 
the  number  of  times  that  a  zero  is  2  spaces  from  another  zero,  etc. 

The  "textural  features"  that  can  be  determined  from  these  GLCM  are  at 
least  28  in  number;  however,  there  are  only  4-6  that  have  been  foimd  to  be 
most  useful  (as  stated  above).  The  mathematical  formulations  and  the 
labeling  used  in  the  table  headings  found  in  the  main  body  of  the  report,  for 
some  of  these,  are  given  below  (these  are  the  ones  which  we  have  found  to 
be  most  usefiil): 


homogeneity  or  energy  or  angular  second  moment: 


N,  N,  f 


/,  =  sz 

1=1  >1  V 


fjui) 

R 


contrast  or  inertia: 


T  p(i  /) 

#s=l  j=\ 


R 


correlation: 


/3  = 


1=1  7=1 


iji.j) 

R 


<y,<y,. 


entropy: 

R 

and  the  definitions  for  the  variables  in  the  above  equations  are: 

P{i,j)  =  the  value  of  the  i,j  element  of  the  GLCM, 

R  s  normalization  constant  (equal  to  the  number  of  neighboring 
resolution  cell  pairs  used  in  computing  the  matrix), 

Ng  =  number  of  gray  levels, 

/=1  7=1 

1=1  i=l 


1=1  7=1 

j=\  i-\ 
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Appendix  D 


Histograms  of  Texture  Metrics 


Figure  D-1.  One  of  the  histograms  used  to  derive  the  parameters  listed  in 
Table  3  of  the  main  body  of  the  report.  This  particular  histogram  shows  the 
distribution  of  the  metric  “correlation  length”  as  applied  to  the  texture  map 
generated  by  the  Stoksik  (1995)  technique. 


Figure  D-2.  One  of  the  histograms  used  to  derive  the  parameters  listed  in 
Table  3  of  the  main  body  of  the  report.  This  particular  histogram  shows  the 
distribution  of  the  metric  “correlation  length”  as  applied  to  the  texture  map 
generated  by  the  Lohmann  (1995)  technique. 


Metric  Value 


Figure  D-3.  One  of  the  histograms  used  to  derive  the  parameters  listed  in 
Table  3  of  the  main  body  of  the  report  This  particular  histogram  shows  the 
distribution  of  the  metric  “correlation  length”  as  applied  to  the  texture  map 
generated  by  the  Midpoint  Replacement  technique. 


Metric  Value 


Figure  D-4.  One  of  the  histograms  used  to  derive  the  parameters  listed  in 
Table  3  of  the  main  body  of  the  report.  This  particular  histogram  shows  the 
distribution  of  the  metric  “correlation  length”  as  applied  to  the  texture  map 
generated  by  the  IRMA  (smoothed  output)  technique. 
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Figure  D-7.  One  of  the  histograms  used  to  derive  the  parameters  listed  in 
Table  3  of  the  main  body  of  the  report.  This  particular  histogram  shows  the 
distribution  of  the  metric  “GLCM  Homogeneity”  as  applied  to  the  texture  map 
generated  by  the  Lohmann  (1995)  technique. 


Figure  D-8.  One  of  the  histograms  used  to  derive  the  parameters  listed  in 
Table  3  of  the  main  body  of  the  report.  This  particular  histogram  shows  the 
distribution  of  the  metric  “GLCM  Homogeneity”  as  applied  to  the  texture  map 
generated  by  the  Midpoint  Replaeement  technique. 
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Figure  D-9.  One  of  the  histograms  used  to  derive  the  parameters  listed  in 
Table  3  of  the  main  body  of  the  report.  This  particular  histogram  shows  the 
distribution  of  the  metric  “GLCM  Homogeneity”  as  applied  to  the  texture  map 
generated  by  the  IRMA  (Smoothed  Output)  technique. 


Figure  D-10.  One  of  the  histograms  used  to  derive  the  parameters  listed  in 
Table  3  of  the  main  body  of  the  report.  This  particular  histogram  shows  the 
distribution  of  the  metric  “GLCM  Homogeneity”  as  applied  to  the  texture  map 
generated  by  the  IRMA  (Unsmoothed  Output)  technique. 
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Figure  D-11.  One  of  the  histograms  used  to  derive  the  parameters  listed  in 
Table  3  of  the  main  body  of  the  report.  This  particular  histogram  shows  the 
distribution  of  the  metric  “GLCM  Correlation”  as  applied  to  the  texture  map 
generated  by  the  Stoksik  (1995)  technique. 


Figure  D-12.  One  of  the  histograms  used  to  derive  the  parameters  listed  in 
Table  3  of  the  main  body  of  the  report.  This  particular  histogram  shows  the 
distribution  of  the  metric  “GLCM  Correlation”  as  applied  to  the  texture  map 
generated  by  the  Lohmann  (1995)  technique. 
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Metric  Value  Histogram 
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Metric  Value 


Figure  D-13.  One  of  the  histograms  used  to  derive  the  parameters  listed  in 
Table  3  of  the  main  body  of  the  report.  This  particular  histogram  shows  the 
distribution  of  the  metric  “GLCM  Correlation”  as  applied  to  the  texture  map 
generated  by  the  Midpoint  Replacement  technique. 


Figure  D-14.  One  of  the  histograms  used  to  derive  the  parameters  listed  in 
Table  3  of  the  main  body  of  the  report.  This  particular  histogram  shows  the 
distribution  of  the  metric  “GLCM  Correlation”  as  applied  to  the  texture  map 
generated  by  the  IRMA  (Smoothed  Output)  technique. 
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Figure  D-15.  One  of  the  histograms  used  to  derive  the  parameters  listed  in 
Table  3  of  the  main  body  of  the  report.  This  particular  histogram  shows  the 
distribution  of  the  metric  “GLCM  Correlation”  as  applied  to  the  texture  map 
generated  by  the  IRMA  (Unsmoothed  Output)  technique. 


Figure  D-16.  One  of  the  histograms  used  to  derive  the  parameters  listed  in 
Table  3  of  the  main  body  of  the  report.  This  particular  histogram  shows  the 
distribution  of  the  metric  “GLCM  Entropy”  as  applied  to  the  texture  map 
generated  by  the  Stoksik  (1995)  technique. 
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Metric  Value  Histogram 


Figure  D-17.  One  of  the  histograms  used  to  derive  the  parameters  listed  in 
Table  3  of  the  main  body  of  the  report.  This  particular  histogram  shows  the 
distribution  of  the  metric  “GLCM  Entropy”  as  applied  to  the  texture  map 
generated  by  the  Lohmann  (1995)  technique. 


Figure  D-18.  One  of  the  histograms  used  to  derive  the  parameters  listed  in 
Table  3  of  the  main  body  of  the  report.  This  particular  histogram  shows  the 
distribution  of  the  metric  “GLCM  Entropy”  as  applied  to  the  texture  map 
generated  by  the  Midpoint  Replacement  technique. 
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Figure  D-19.  One  of  the  histograms  used  to  derive  the  parameters  listed  in 
Table  3  of  the  main  body  of  the  report.  This  particular  histogram  shows  the 
distribution  of  the  metric  “GLCM  Entropy”  as  applied  to  the  texture  map 
generated  by  the  IRMA  (Smoothed  Output)  technique. 


Melric  Value 


Figure  D-20.  One  of  the  histograms  used  to  derive  the  parameters  listed  in 
Table  3  of  the  main  body  of  the  report.  This  particular  histogram  shows  the 
distribution  of  the  metric  “GLCM  Entropy”  as  applied  to  the  texture  map 
generated  by  the  IRMA  (Unsmoothed  Output)  technique. 
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Figure  D-21.  One  of  the  histograms  used  to  derive  the  parameters  listed  in 
Table  3  of  the  main  body  of  the  report.  This  particular  histogram  shows  the 
distribution  of  the  metric  “Fractal  Dimension”  as  applied  to  the  texture  map 
generated  by  the  Stoksik  (1995)  technique. 
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Figure  D-22.  One  of  the  histograms  used  to  derive  the  parameters  listed  in 
Table  3  of  the  main  body  of  the  report.  This  particular  histogram  shows  the 
distribution  of  the  metric  “Fractal  Dimension”  as  applied  to  the  texture  map 
generated  by  the  Lohmann  (1995)  technique. 


Metric  Value 


Figure  D-23.  One  of  the  histograms  used  to  derive  the  parameters  listed  in 
Table  3  of  the  main  body  of  the  report.  This  particular  histogram  shows  the 
distribution  of  the  metric  ‘‘Fractal  Dimension”  as  applied  to  the  texture  map 
generated  by  the  Midpoint  Replacement  technique. 


Metric  Value 


Figure  D-24.  One  of  the  histograms  used  to  derive  the  parameters  listed  in 
Table  3  of  the  main  body  of  the  report.  This  particular  histogram  shows  the 
distribution  of  the  metric  “Fractal  Dimension”  as  applied  to  the  texture  map 
generated  by  the  IRMA  (Smoothed  Output)  technique. 
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Figure  D-25.  One  of  the  histograms  used  to  derive  the  parameters  listed  in 
Table  3  of  the  main  body  of  the  report.  This  particular  histogram  shows  the 
distribution  of  the  metric  “Fractal  Dimension”  as  applied  to  the  texture  map 
generated  by  the  IRMA  (Unsmooothed  Output)  technique. 
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Figure  D-26.  One  of  the  histograms  used  to  derive  the  parameters  listed  in 
Table  3  of  the  main  body  of  the  report.  This  particular  histogram  shows  the 
distribution  of  the  metric  “GLCM  Contrast”  as  applied  to  the  texture  map 
generated  by  the  Stoksik  (1995)  technique. 


Figure  D-27.  One  of  the  histograms  used  to  derive  the  parameters  listed  in 
Table  3  of  the  main  body  of  the  report.  This  particular  histogram  shows  the 
distribution  of  the  metric  “GLCM  Contrast”  as  applied  to  the  texture  map 
generated  by  the  Lohmann  (1995)  technique. 


Figure  D-28.  One  of  the  histograms  used  to  derive  the  parameters  listed  in 
Table  3  of  the  main  body  of  the  report.  This  particular  histogram  shows  the 
distribution  of  the  metric  “GLCM  Contrast”  as  applied  to  the  texture  map 
generated  by  the  Midpoint  Replacement  technique. 
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Figure  D-29.  One  of  the  histograms  used  to  derive  the  parameters  listed  in 
Table  3  of  the  main  body  of  the  report.  This  particular  histogram  shows  the 
distribution  of  the  metric  “GLCM  Contrast”  as  applied  to  the  texture  map 
generated  by  the  IRMA  (Smoothed  Output)  technique. 


Figure  D-30.  One  of  the  histograms  used  to  derive  the  parameters  listed  in 
Table  3  of  the  main  body  of  the  report.  This  particular  histogram  shows  the 
distribution  of  the  metric  “GLCM  Contrast”  as  applied  to  the  texture  map 
generated  by  the  IRMA  (Unsmoothed  Output)  technique. 
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Appendix  E 

AR  (1)  Process:  A  Brief  Discussion 


According  to  Kendall  (1976),  an  autoregressive  process  is  a  stationary 
series  which  is  “...generated  by  a  mechanism  in  which  the  value  of  the 
series  at  time  t  is  expressed  in  terms  of  the  past  values  ~  a  systematic 
dependence  on  past  history  -  plus  a  ‘disturbance’  term  f  happening  at  time 
r.”  Stationarity  is  defined  here  as  applying  to  a  series  which  has  had  its 
trend  removed  or  has  never  had  a  trend.  The  autoregressive  process  of 
order  p  is  one  defined  by  : 

u,  =  - ^p^i-p  +  • 

Because  the  form  of  this  equation  suggests  a  regression  process,  the  term 
“autoregressive”  has  been  coined  -  in  fact,  this  is  not  necessarily  so.  The 
first  order  autoregressive  process,  AR(1),  is  also  know  as  a  Markov  process 
and  is  defined  by: 

u,  +f,. 

Rewriting  this  (still  following  Kendall)  as: 

U,  =PU,_,  +€, 

then  allows,  after  rearranging  of  terms  and  further  calculations,  for  the 
determination  that  the  autocorrelation  of  u,  at  lag  1  (p)  is  the  same  as  the 
value  of  the  coefficient  -a^  in  the  above  equation.  The  limits  on  the  value 
of  p  are  that  it  be  less  than  or  equal  to  1  and,  if  the  autocorrelation  function 
is  to  decay  without  oscillation,  then  it  should  also  be  positive;  otherwise,  the 
autocorrelation  function  will  oscillate  as  it  decays. 

Another  way  of  considering  the  AR(1)  process  is  that  it  is  a  linear  process 
which  has  as  its  input  white  noise  (Jenkins  and  Watts,  1969)  which  is 
filtered  through  the  equivalent  of  an  RC  circuit  whose  transfer  fimction  is 
given  by: 

h{v)  =  je 
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Such  a  linear  process  is  given  by: 


t^+{x(i)-m)  =  z{') 

where: 

X{t)  =  system  output 
Z{t)  =  system  input  (white  noise) 

T  =  system  time  constant 
fd  -  mean  of  the  process. 

This  has  a  digital  analog  which  is  just  that  described  above  and  given  by 
Kendall. 
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Appendix  F 


Mid-Point  Displacement  Code 
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fractscn.c  implements  the  mid-point  displacement  pseudocode  found  on  page  100  of  the  book  "The  Science  of  Fractal 
Images"  by  Peitgen  et  al. 
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FILE  *outptrl; 
float  huge  *  Y ; 


int  maxlevel;//maximal  #  of  recursions, N=2^maxlevel 
float  sigma;//initial  standard  deviation 
float  H;//parameter  H  determines  fractal  dimension  D=3-H 
int  addition://parameter(0/!=0)tums  random  additions  offron 
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printf("\n"); 

printf("  Size  of  image  =  N  by  N  where  N=2^{maxlevel}\n"); 


if(  (outptrl=fopen("mpscn.pix","w+b"))  =  NULL) 
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/*call  function  to  generate  fractal  image  */ 
MidPointFM2D(Y,maxlevel,sigma,H, addition, seed); 


/*calculate  max,  min,  mean  and  variance  of  N  by  N  subimage*/ 
sum=0; 

for(i=0;i<N;++i){ 

for(j=0;j<N;++j){ 
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/*output  square  N  by  N  array*/ 
for(i=0;i<N;++i){ 


for(j=0;j<N;++j){ 

Y[i*(N+l)+j]  =  (Y[i*(N+l)+j]-tmean)*stdev/tstdev+mean; 
if(Y[i*(N+l)+j]  <  0.)Y[i*(N+l)+j]=0.; 

if(Y[i*(N+l)+j]  >  (float)max_pix_val  )Y[i*(N+l)+j]=(float)maxjpix_val; 
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for(i=0;i<N;++i){ 

for(j=0;j<N;++j){ 

temp=(double)(Y[i*(N+ 1  )+j]-tmean); 


tvar=(float)(sumsq/(  (double)N*(double)N-l. )); 
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int  i,N,stage; 

float  delta;//real  holding  standard  deviation 
int  x,y,yO,D,d;//integer  array  indexing  variables 


srand(seed); 

N=pow(2,maxlevel); 

/*set  the  initial  random  comers*/ 
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/*displace  other  points  also  if  needed*/ 
if(addition){ 

for(x=0;x<=N;x+=D)  { 


for(y=0;y<=N;y+=D){ 

X[x*(N+l)+y]  +=  delta*  gaussian(0.,l.); 
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for(y=d;y  <=N-d;y+=D)  { 

X[x*(N+l)+y]  =  f4(delta,X[x*(N+l)+y+d],X[x*(N+l)+y-d],\ 

X[(x+d)*(N+l)+y],X[(x-d)*(N+l)+y]); 


}/*end  X  for*/ 

/*displace  other  points  also  if  needed*/ 
if(addition){ 
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//Initialize  rand  by  srand  before  this  function  is  called. 


double  pi=3.14159265358979; 

u=(float)((float)rand()/(float)RAND_MAX); 

if(u=0.)u=.000000001; 
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Appendix  G 


GLCM  Code 


glcmtex.c  implements  the  algorithm  described  in  the  Gabriele  Lohmann  article 
"Analysis  and  Synthesis  of  Textures;  A  Co-Occurrence-Based  Approach". 
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void  find_hist(int  huge  *,int  huge  *,int  huge  *,int  huge  *,int  huge  *,int); 
void  ran_image(int  huge  *,int  huge  *,int,int,int); 
long  int  metric(int  huge  *,int  huge  *,int); 


long  int  change(int  huge  *,int  huge  *,int  huge  *,int,int,int,int,int,int,int,int,int); 
void  update_glcm(int  huge  *,int  huge  *,int,int,int,int,int,int,int,int,int); 
void  out_mess(char  [],FILE  *, int, int, int  huge  *); 
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/*  the  number  of  rows  and  columns  are  fixed  at  32.  */ 


/♦greylevels  set  to  128  to  limit  the  size  of  the  huge  arrays  h0,h45,etc. 
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hl35=(int  huge  *)farcalloc(16384UL,(unsignecl  long)sizeof(int)); 
tO=(inthuge  *)farcalloc(l  63  84UL, (unsigned  long)sizeof(int)); 
t45=(inthuge  *)farcalloc(16384UL,(unsigned  long)sizeof(int)); 
t90=(inthuee  *)farcalloc(l 63  84UL, (unsigned  long)sizeof(int)); 
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farfree((int  far  *)t0); 
farfree((int  far  *)t45); 
farfree((int  far  *)t90); 


farfree((int  far  *)tl35); 
farfree((int  far  *)hist); 
farfree((int  far  *)texture); 
fclose(infptrl); 


C/3 

"o 

JJ 

T3 

oj 

s 

N 

TO 

‘So 

C/3 

.S 

1 

1/3 

1/3 

o 

o 

feb 

c/3 

cd 

X3 

X 

x: 

o 

1/3 

X 

13 

> 

> 

(D 

bD 

(U 

cb 

C 

H 

c/3  O 

^  B 
*>»  9 

u  ♦  'O 

feb  a>  * 

^  00  bp 

c  S 

-a  S 

-fi 

s  .s  ^ 
^  :5  ^ 

SP  «>  So 

3  00  3 

B  £ 

•B 


S^r>  "/-T 
op  O  fD 
H  On  ^ 

S 

r  *  * 

flj  O  ^ 

^  bO  Gp 

.1 4  -1 

(yf  irT  O 
^  ^  On 

o  5p  * 

^  f 
s  s  .s 


o 

X 

V  J 

§ 

* 

* 

c 

’g 

Sti. 

‘S 

bO 

bO 

<U 

cd 

:b 

:b 

c/3 

g 

U) 

X 

X 

o 

Jo 

no 

B 

.S 

Oh 

c 

o 

* 

rg 

'o 

> 

Oh 

Oh 

X 

O 

o 

bo 

CO 

O 

cd 

X 

o 

•  imH 

9--r 

* 

Oh 

'S 

-S 

O 

e 

c/3 

CO 

nr? 

.5  O 

’O  i-J 

*§  <o 
■g  eS 

sS 

a  ON 


^  o 

=  m 

"S 

O 
cd  ON 


S  •' 

-2  o 

■g  M 


=5  'C 

.  §  S 

O  O 

S  -cj  -c 

I  H  H 


(U 

tS 

^  bo 

^  c 


.S2  o 

x:  T3 


122 


intrs0[256],rs45[256],rs90[256],rsl35[256]; 

intbr[256]; 

int  old_delta0,old_delta45  ,old_delta90,old_delta  135; 
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/*find  initial  random  image  "texture"  with  thb  histogram*/ 


srand(seed); 

ran_image(hist,texture,greylevels, rows, cols); 
out_mess("random  texture  start", outptrl, rows, cols, texture); 
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for(i=0;i<rows*cols;++i){ 

/*pick  two  distinct  pixels  to  swap*/ 
pick2(rows,cols,&i  1  ,&j  1  ,&i2,&j2); 
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delta  increases  as  a  result*/ 

/*call  to  update  glcm  must  proceed  swap  of  pixels*/ 
update_glcm(texture,tO,greylevels,rows,cols,i  1  ,j  1  ,i2,j2,0, 1 ); 


update_glcm(texture,t45,greylevels, rows, cols, i  1 J 1 ,12  j2, 1,1); 
update_glcm(texture,t90,greylevels,rows,cols,i  1 J 1 ,12  J2, 1 ,0); 
update_glcm(texture,t  1 3  5 , grey  levels, rows, cols, 1 1 J 1 ,12  J2, 1,-1); 
swap(texture,cols,rows,l  1 J 1 ,12  J2); 
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T  *=  .99; 

}whlle(enough==0);  /*end  do  while  loop*/ 
}/*endtxtglcml  function*/ 


void  out_mess(char  message[80],FILE  *outptr,int  rows,int  cols,int  huge  *target){ 
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printf("%s\n", message); 
for(i=0;i<rows;++i){ 


for(j=0;j<cols;++j){ 

printf("%d  ",target[i*rows+j]); 
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char  * 


for(i=0;i<rows*cols;++i){ 

fscanf(inQ)tr,"%s",v); 

target[i]=atoi(v); 
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int  ij,gl,glplus,jstart,jend; 
for(i=0;i<grey  levels*  greylevels;++i)  { 
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void  rowsums(int  huge  *  input, int  huge  *  output, int  grey  levels)  { 

/*  input  is  a  2d  square  matrix  of  size  grey  levels*  grey  levels 

output  is  a  vector  of  size  greylevels  containing  the  row  sums  of  input 


for(i=0;i<greylevels;++i){ 

outputril=0; 
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void  ran_image(int  huge  *hist,int  huge  *texture,int  grey  levels, int  rows,int  cols){ 


/*  This  function  generates  a  random  rows  by  cols  image  with  grey  level  histogram  "hist". 
The  random  number  initialization  fiinction  srand()  has  already  been  called 
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long  int  delta; 
delta==^0; 


for(i==0;i<grey  levels*  grey  levels;++i)  { 
delta  +=  (int  long)abs(ml[i]-m2[i]); 
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/*  This  function  calculates  the  value  "delta"  resulting  from  a  swap  of 
the  greylevel  values  at  pixels  01, if)  and  02,i2).  "jl"  is  a  row  index. 


Delta"  is  the  change  in  the  "distance"  between  the  target  and  texture  GLCMs 
before  and  after  the  swap.  So  if  "delta"  is  positive,  the  swap  has  moved  the 
texture  GLCM  farther  away  from  the  target  GLCM. 
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Pm=j  1  m*  cols+i  1  m; 

Pp=jlp*cols+ilp; 

Q=j2*cols+i2; 


Qm=j2m*cols+i2m; 

Qp=j2p*cols+i2p; 

/*find  the  values  occurring  in  the  target  matrix*/ 
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if(in_bounds(i2m,j2m,cols,rows)){ 

check_add(address,deltat,&index,(long)(d*g+a),- 1 ); 
check_add(address,deltat,&index,(long)(d*  g+b),  1 ); 


else  if(i2=ilp&&j2=jlp){ 

if(in_bounds(i  1  m,j  1  m, cols, rows))  { 
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check_add(address,deltat,&index,(long)(b*  g+c),- 1 ); 
check_add(address,deltat,&index,(long)(e*  g+c),  1 ); 


i_bounds(i2m,j2m,cols,rows))  { 
check_add(address,deltat,&index,(long)(d*g+e),- 1 ); 
check_add(address,deltat,&index,(long)(d*  g+b),  1 ); 
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int  i  1  m,i  1  p,i2m,i2p  J 1  mj  1  p,j2m  J2p,P,Pm,Pp,Q,Qm,Qp,a,b,c,d,e,f; 


int  long  delta; 
ilm=il-dx; 
jlm=jl-dy; 
ilp=il+dx; 
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/♦cover  all  the  three  cases*/ 
if(i2=i  1  m&&j2=j  1  m){ 

if(in_bounds(i  1  m,j  1  m, cols, rows))  { 


if(in_bounds(i  1  p,j  1  p, cols, rows))  { 
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-t[e*g+f]; 

++t[d*g+fl; 


if(in_bounds(i  1  m  j  1  m, cols, rows))  { 
-t[a*g+b]; 
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/*  This  ftinction  determines  whether  the  pixel  with  calculated  coordinates 
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/*  This  fiinction  picks  two  pixels  at  random  */ 


il=rand()  %  cols; 
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br(i=0;i<*index;++i){ 

if(address[i]=address_item)  { 
deltat[i]  +=  change; 


return; 
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Appendix  H 


Correlation  Length  Code  (IRMA) 


This  C  code  implements  the  SWOE  texture  generation  eode  based  on  the  IRMA  model. 

This  code  is  a  translation  of  the  Pascal  code  in  SWOE  Report  91-3,  "Generating  Textures  for  Synthetic  Thermal 
Seenes",  by  B.M.  Sabol  and  L.K.  Baliek. 


(D  ,  I 
^  fa 

/I  \  ^ 


u  r- 

(N 


(U 

s 

oo 

> 

m 

<u 

c/5 

cn 

C/5 

o^ 

"o 

o 

'O 

c/5 

2 

o 

VO 

Oh 

(L> 

o 

c 

m 

in 

'O 

00 

c 

T— H 

cd 

6 

fN 

g 

m 

o 

o 

CN 

too 

m 

C 

1 

>< 

0) 

JO 

>< 

s 

<u 

c 

C 

H 

AAA 

j3  j=:  -a 

^  .2  -S 
—  ^  2 
«  -s  S 

V  V  V 


^ 


A 

£  ■§ 
’S  s 


/ — ^  * 

*  <D 
8)  ^ 

^  I-* 

*  .S 

OD  'S  .5 
p  .a 

-S 

C  •'  * 
'-S'  pj 

S  ^  3 
•S 

W  ja  ctf 

E  o  w 
w  w  w 

(U  GO  W 

M  ??  2 


V 

<U 

V 

<U 

V 

IE 

•^1 

p 

T3 

o 

cd 

o 

c/5 

"o 

73 

'g 

rs 

td 

g 

g 

’o 

'o 

*o 

O 

> 

> 

> 

pH 

146 


float  gasdevO; 

void  autoregression(int  huge  *,int,int, float, float, float, float); 
void  out_char(FILE  *,int,int,int  huge  *); 


iving_average(int  huge  *,int,int); 
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h0=(int  huge  *)farcalloc(16384lJL,(unsigned  long)sizeof(int)); 
texture=(int  huge  *)farcaIloc(16384lJL,(unsigned  long)sizeof(int)); 
printf("\nMemory  left  on  the  far  heap:  %lu\n",farcoreleft()); 


rows=atoi(argy  [  1  ]); 
cols=atoi(arevr21); 
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fclose(outptrl); 

fclose(outptr2); 

fclose(outptr3); 


}/*end  main*/ 
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hold  1  [lpv]=one_ninth; 
one_ninth=0.0; 


texture[(lph-l)*cols+tt]=hold2[tt]; 
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void  autoregression(int  huge  *texture,int  rows,int  cols,\ 


int  ij,h; 

float  rij,rijeq,zerozero,adis  1  ,adis2,rmean,tl  ,t2,rstd; 
float  pixelprev,pixelabov,pixelupbak,adjust; 
float  parrr2ir442]; 
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pixelabov=parr[0]|j  ]; 

pixelupbak=parr[0][j- 1  ]; 

adjust=adis  1  *pixelprev+adis2*pixelabov\ 


-adis  1  *adis2*pixelupbak+rijeq; 
parr[l]|j]=adjust; 
if(adjust<(float)  1  .)adjust=  1 
if(adjust>(float)max_pix_val)adjust=max_pix_val; 
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glgset=vl*fac; 

gasdev=v2*fac; 

gliset=l; 


gliset=0; 
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glix  1  =(ic  1  -idum)%m  1 ; 
glix  1  =(ia  1  *  glix  1 +ic  1  )%m  1 ; 
glix2=glixl  %  m2; 
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printf("%d  ",target[i*cols+j]); 


void  in_image(FILE  *infptr,int  rows,int  cols,int  huge  *target){ 


* 


a 


u 

c/3 

< 

C 

X 

B 

S 

a 

B 

o 

C/3 

"3 

Oh 

.S 

c 

V-> 

o 

c 

a 


o 

o 

* 

I/) 

O 

>  o 
*  II 

.S  'o  ^ 


(/3  . 

vp 

^  > 
X  C& 

±  .S  rtU 

s  ^ 

s  g) 


x: 

H 

* 


bX) 


* 


o 

o 


c/3 

I 

fM 


d 

a 

t 

o 

* 


3 

'o 

> 


* 

s 

cd 

O 

'4-» 

■4-* 

o 

§ 

hO 

o 


(U 

§) 

<+H 

O 

X 


cd 

s 

cd 


d 

& 

§ 

d 

o 

d 

a 


c/3 


o 

o 

o 

o 


d  .. 
'6  4 


O 

O 

o 

o 


if 


3  .- 
3  3 


bO  fa 

s  ^ 

^  5  I 

cd  "d  "d 
C  (U  <u 

.| 

l/J  c/3 

d  d 
.S  3  3 


156 


for(i=0;i<rows;++i){  /*find  max  and  min  pixel  value*/ 
for(j=0y<cols;-H-j){ 


numpix=i*cols+j ; 

if(target[numpix]>maxt)maxt=target[numpix]; 

if(target[numpix]<mint)mint=target[numpix]; 
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Appendix  I 


“Accurate  Fractal”  Code 


/*July  5,1995,accfrac.c:This  program  implements  the  algorithm  found  in 
the  article  "Practical  Synthesis  of  Accurate  Fractal  Images"  by 
Stoksik  et  al.  in  Graphical  Models  and  Image  Processing  Vol  57, 

No.  3,pp206-219,1995. 

It  includes  the  fimctions  AccFB2d  and  Gaussian 
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float  huge  *Y; 

int  maxlevel;//maximal  #  of  recursions, N=2^maxlevel 

float  D_image;//width  of  the  image 

float  H;//parameter  H  determines  fractal  dimension  D=3-H 


const  int  max_pix_val=127; 

unsigned  int  seed;//seed  value  for  random  number  generator 

unsigned  char  ucharxij; 

float  maxxij=-l  000000., minxij=l  000000.; 
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maxlevel=atoi(argv[  1  ]); 

N=pow(2,maxlevel); 

printf("Memory  in  far  heap:  %lu\n",farcoreleft()); 

Y=(float  huge  *)farcalloc(((long)N+l)*((long)N+l),sizeof(float)); 
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fclose(outptrl); 
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/*interpolate  and  offset  interior  grid  points*/ 
sigma_eps  *=  pow(.5,.5*H); 
for(x=d;x<=N-d;x+=D)  { 

for(y=D;y<=N-d;y+=D)  { 

X[x*(N+l)+y] 
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float  gaussian(float  mean, float  stdev) 

//Before  this  function  is  called  rand  needs  to  be  initialized  by  srand 


float  u,v,x; 

double  pi=3. 14159265358979; 

u=(float)((float)rand()/(float)RAND_MAX); 

if(u==0.)u=.000000001; 


^  11 
^  X 

|j 

c/3 

^  o 

> 

td 

X 

O  11 

4h' 

o 

c 

03  II 

U 

II 

w 

II 

B 

w 
li  ^23^ 

>  Jt3 

g 

1!< 

u 

CIh 


c/3 

m  w 

i  s 

cd  W) 

13 

"Td 

•5  i 

o  5 
•'  m 

I? 

o' 

X  + 
^  o 

cd  X 
O 
03 
cd 

4-*  ^ 

o 

03 

'S 

o 

03 


e 


a 

cd 

c/3 

W) 

* 

CT  T3 

X  5 
• 

o 

'n  ^ 
o  ^ 
X  + 

^4mJ 

cd  X 
o 

'n  g 

S  I 

U 

o 

q:3 


a 

'S 

o 

03 


166 


Appendix  J 


Discussion  of  Fractal  Dimensions 


There  are  several  points  which  must  be  considered  under  the  topic  of 
“fractal  analysis”  and/or  synthetic  scene  generation  using  a  “fractal  seed”. 
These  can  be  sununarized  to  three  main  categories: 

1 .  Underlying  Model:  Is  the  texture  monofractal  or  multifractal? 

2.  Metrics:  Which  tool  does  one  use  to  measure  the  “dimension”  (e.g., 
power  spectrum  or  variogram)? 

3.  Generation/Simulation  Mechanism:  Which  synthetic  “scene”  tool  to  use 
(e.g.,  midpoint  displacement  or  “accurate  fractal”)? 

For  the  purposes  of  the  study  reported  here,  we  assumed  a  monofractal 
model  and  used  the  power  spectrum  as  the  tool  for  determining  the  “fractal 
dimension”  of  the  texture.  We  then  chose  the  midpoint  displacement 
algorithm  as  the  fractal  “scene”  generation  tool.  As  part  of  the  initial  study, 
we  also  applied  another  scene  generation  technique  (Stoksik,  1995)  as  it 
was  suggested  that  it  could  more  accurately  replicate  a  texture  whose  fractal 
dimension  were  the  same  as  the  “input”  dimension  —  this  part  of  the  effort 
was  undertaken  after  the  main  study  was  already  begun;  hence,  it  was 
included  only  in  appendices  to  the  main  report  (Appendix  A:  Modified 
Mid-Point  Displacement  Algorithm)  and  even  then,  not  in  as  complete  a 
manner  as  the  other  techniques.  The  approach  taken  in  this  report  has  left 
some  unanswered  questions  which  it  was  felt  should  be  addressed  even  if 
they  could  not  be  answered.  What  instigated  this  discussion  is  that  the 
“measured”  fractal  dimensions  of  some  of  the  texture  maps  were  outside  the 
accepted  bounds  for  that  metric;  i.e.,  fractal  dimensions  greater  that  3  were 
obtained  for  some  surfaces.  Post  facto  efforts  to  rationalize  these  results 
have  arrived  at  some  conclusions  which  are  of  interest.  These  will  be 
addressed  below,  briefly. 

First,  there  are  techniques  to  determine  whether  or  not  the  object  of  study  is 
a  fractal  --  we  chose  to  not  use  them  at  the  outset  for  several  reasons: 

1.  limited  resources, 

2.  an  incomplete  understanding,  on  our  part,  of  this  whole  fractal  business, 
and 

3.  more  than  likely,  there  are  no  accepted  procedures  to  follow  anyway 
(Tsonis  and  Eisner,  1995) 


Some  of  these  are  addressed  in  a  paper  by  Lovejoy  et  al.  (1995).  For 
example,  common  techniques  to  measure  the  fractal  dimension  of  a 
“process”  use  variogram  or  spectral  methods  to  arrive  at  a  “fractal 
dimension”  -  according  to  Lovejoy  et  al.(1995),  these  methods  “...measure 
the  scaling  exponent  of  the  second  order  moments... which  is  not  a  fractal 
dimension”  unless  monofractality  is  assumed  and,  indeed,  is  true. 
Otherwise,  the  measurement  is  erroneous.  Many  making  “fractal 
measurements”  discuss  complex  techniques;  but  then,  choose  those  ways 
of  doing  things  with  which  they  are  most  comfortable  and  which  are  in 
many  cases  the  simplest  (we  were  no  different  in  that  we  chose  techniques 
that  were  already  in  use  by  us  and  with  which  we  were  already 
comfortable).  These  simple  techniques  may  work  for  truly  fractal 
phenomena;  however,  when  the  process  is  multifractal,  much  more 
complicated  techniques  are  required.  Those  who  shun  power  spectral 
methods  because  of  their  complexity  will  find  multifractal  analysis  much 
more  difficult  to  both  accomplish  and  understand.  “A  fractal  is  a 
geometrical  set  of  points;  a  multifractal  is  a  mathematical  measure.  ...many 
results  derived  for  fractal  sets  and  monofractal  functions  will  not  apply 
(Lavallee  et  al.,  1993).”  So,  were  the  effort  described  in  this  report  to  be 
done  again,  it  is  highly  recommended  that  this  issue  between  fractal  and 
multifractal,  and  even  between  fractal  and  Euclidean,  be  resolved  first.  The 
point  here  is  that  if  our  target  texture  were  non-fractal  or  even  multifractal, 
then  the  results  of  applying  our  metrics  have  no  meaning  and  the  validity  of 
the  number  used  by  the  synthetic  scene  generation  algorithm  would  then  be 
suspect. 

Second  is  the  problem  of  the  metrics.  It  turns  out  that  there  are  a  large 
number  of  metrics  available  to  determine  the  fractal  dimension  of  an 
“object”  of  study  --  some  are  more  or  less  appropriate  depending  on  the 
object  while  others  are  not  and  still  others  are  misunderstood  (e.g.,  see 
Kinsner,  1994).  For  example,  there  is  the  subject  of  using  the  power 
spectrum  as  a  tool  to  determine  the  fractal  dimension  of  a  surface  —  Voss 
(1988)  gives  equations  for  and  discusses  the  use  of  the  1-,  2-,  and  3- 
dimensional  power  spectrum  for  the  determination  of  the  fractal  dimension 
of  those  “objects”:  a  line,  a  surface,  and  a  volume,  respectively.  He 
describes  the  relationship  between  them  with  the  conclusion  on  many 
analysts’  part  that  one  can  determine  the  power  spectrum  of  a  vertical  cut 
through  a  surface  and  that  the  slope  of  this  is  just  one  less  than  the  power 
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spectrum  of  the  surface  (2-dimensional)  —  it  may  be  that  this  conclusion  is 
erroneous  (Clarke  and  Schweizer,  1991))  and  that  the  correct  relationship 
between  the  fractal  dimension  of  a  surface  and  a  cut  through  that  surface  is 
for  a  horizontal  cut  --  not  a  vertical  cut!  This  then  begs  the  question  of  how 
one  determines  the  power  spectrum  of  a  line  which  forms  the  horizontal  cut 
or  is  it  truly  acceptable  to  compare  results  from  vertical  cuts  with  2- 
dimensional  power  spectra?  Huang  and  Turcotte  (1990)  and  Turcotte 
(1992)  seem  to  arrive  at  acceptable  results  this  way  -  they  describe  the 
results  of  applying  power  spectra  to  small  portions  of  the  topographic  map 
of  Oregon  with  good  results  in  the  comparison  between  one-dimensional 
vertical  cut  and  two-dimensional  spectra  (though,  in  an  earlier  paper 
(Huang  and  Turcotte,  1989)  which  studied  the  state  of  Arizona,  the 
agreement  was  not  very  good  and  they  do  not  discuss  the  discrepancy, 
anywhere!).  Lam  and  De  Cola  (1993)  describe  an  isarithm  technique  for 
dealing  with  the  horizontal  cut  as  well  as  a  variogram  technique  for  dealing 
with  the  whole  surface.  We  have  applied  the  variogram  technique  to 
several  of  our  texture  maps  for  the  purpose  of  measuring  with  another 
technique.  This  is  discussed  next. 

The  variogram  technique  which  we  applied  to  11  each  midpoint 
displacement  and  accurate  fractal  maps  does  not  seem  to  clear  up  the 
measurement  problem  and,  in  fact,  may  not  even  be  appropriate  (Lavallee  et 
al.,  1993)  Regardless,  figures  1  and  2  show  examples  of  the  variogram  for 
two  different  texture  maps,  one  each  for  the  two  fractal  generation 
techniques.  The  least  squares  fit  to  the  variogram  is  supposed  to  yield  the 
“Hurst  Parameter”,  H.  For  all  of  the  variograms  used  here,  we  fitted  only 
over  the  “linear”  portion  with  the  following  results  as  given  in  Table  1.  So, 
we  have  a  paradox:  of  the  two  simulation  techniques,  the  one  which  is  the 
most  “accurate”,  performs  poorest.  At  least,  this  is  so  for  the  variogram 
measurement  algorithm;  perhaps,  another  metric  would  yield  different 
results  and,  the  sample  size  may  be  too  small  to  derive  a  valid  conclusion 
on  this  point,  anyway.  According  to  Tsonis  and  Eisner  (1995),  these 
techniques  are  inadequate  to  demonstrate  scaling  —  what  we  have  done  is 
assumed  scaling  and  not  checked  the  viability  of  the  assumption. 


Still  another  issue,  also  discussed  in  the  paper  by  Tsonis  and  Eisner  (1995) 
and  others,  is  that  for  the  type  of  data  with  which  we  are  working,  self- 
affine  processes  are  important.  This  also  means  that  the  dimensions  that  we 
arrive  at  are  supposed  to  be  the  result  of  an  “ensemble  average”  -  the 
sample  that  we  are  working  with  may  or  may  not  be  representative  of  this 
process.  Similarly,  the  synthetic  generation  processes  need  to  be  run  many 
times  to  arrive  at  a  statistically  significant  sample  size  (Tsonis  and  Eisner, 
1995,  in  their  experiment  ran  their  sample  size  to  1000). 


Table  1.  Variogram  measurement  of  the  fractal  dimension  for  11  texture 
maps  for  each  of  the  two  different  “fractal”  texture  generation  algorithms 
(the  algorithms  were  run  with  H=0.28  as  a  seed). 


Source 

Hurst  parameter 
(average) 

Hurst  parameter 
(standard  deviation) 

Fractal  dimension 
P=3-H) 

Midpoint 

0.286 

0.064 

2.714 

displacement 

Accurate  fractal 

0.200 

0.066 

2.800 

As  is  seen  by  the  results  shown  in  Table  1,  different  scene  generation 
techniques  yield  different  results  -  which  one  is  correct?  Or  Best?  Again, 
more  work  is  indicated.  However,  the  problem  which  will  remain  is  one  of 
separating  the  measurement  technique  from  the  process  being  measured. 
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